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ABSTRACT 

We investigate the birth and evolution of Galactic isolated radio pulsars. We begin by estimat- 
ing t heir birth space velocity distribution from proper motion measurements of iBrisken et al.l (120021 
l2003f ). We find no evidence for multimodality of the distribution and favor one in which the absolute 
one-dimensional velocity components are exponentially distributed and with a three-dimensional mean 
velocity of 380lgg km s 1 . We then proceed with a Monte Carlo-based population synthesis, mod¬ 
elling the birth properties of the pulsars, their time evolution, and their detection in the Parkes and 
Swinburne Multibeam surveys. We present a population model that appears generally consistent with 
the observations. Our results suggest that pulsars are born in the spiral arms, with a Galactocentri c 
radial distribution that is well described by the functional form proposed bv lYusifov fe Kiiciikl (|2004f ). 
in which the pulsar surface density peaks at radius ~ 3 kpc. The birth spin period distribution extends 
to several hundred milliseconds, with no evidence of multimodality. Models which assume the radio 
luminosities of pulsars to be independent of the spin periods and period derivatives are inadequate, 
as they lead to the detection of too many old simulated pulsars in our simulations. Dithered radio lu¬ 
minosities proportional to the square root of the spin-down luminosity accommodate the observations 
well and provide a natural mechanism for the pulsars to dim uniformly as they approach the death 
line, avoiding an observed pile-up on the latter. There is no evidence for significant torque decay (due 
to magnetic field decay or otherwise) over the lifetime of the pulsars as radio sources (~ 100 Myr). 
Finally, we estimate the pulsar birthrate and total number of pulsars in the Galaxy. 

Subject headings: pulsars: general — methods: statistical — stars: neutron - stars: kinematics — 
Galaxy: structure 


1. INTRODUCTION 

The birth and evolution of pulsars are of considerable 
interest. The spatial distribution of pulsars at birth may 
be used to associate them with their progenitors. Their 
initial spin periods may also be related to those of the 
progenitors’ cores, which are not well predicted by theory 
due to differential rotation (lEndal & Sofial fl978h . The 
birth properties of neutron stars are also intimately re¬ 
lated to the physics of core-collapse supernovae in which 
most are thought to be formed. For example, their birth 
space velocities, spin periods, and magnetic fields put 
severe constraints on the mechanis ms that may account 
for their high observed velocities ( Laill2003j ). The com¬ 
parison of the Galactic pulsar birth and supernova rates 
may provide further insight into supernovae by quanti¬ 
fying the fraction which leaves behind a neutron star. 
It may also be possible to clarify whether all neutron 
stars unaffected by peculiar conditions (such as accre¬ 
tion from a binary companion or an anomalously high 
magnetic field) are in fact radio pulsars, although not al¬ 
ways beamed toward us. Knowledge of the spin evolution 
of pulsars is particularly valuable in elucidating whether 
magnetic field decay occurs in isolated neutron stars. A 
correlation between the spin down of pulsars and the evo¬ 
lution of their radio luminosity may finally shed light on 
the long-standing problem of the pulsar radio emission 
mechanism. 
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One way to probe the birth and evolution of pulsars 
is to study the population as a whole. While the earlier 
efforts in this field relied on simplified analytical or 
semi-analyti c al treatments fe.g.. iGiinn fc Qs trikei!ll 970l: 
Large! 119711 iTavlor fc Manchester! 119771 : iDavies et al'l 
19771 : ILvne et al.l 1198511 . the trend over the last two 


decades has been to attempt detailed computational 
modelling of the pulsar population and the selection 
effects affecting the observed sample, often including the 
evolution of t he synth etic pulsars from birth to de tec¬ 
tion fe. g.. IStollmanl 1198731: Emmcring & = X'hgygligr 
1980; IBhattacharv a et al.l 1 19921 [ Hartman et al. 


tmattacnarva e t ai.i I luuzt i Hartman etai 
1997 : IL orimer et al.l Il997t I Arzoumanian et al.l 1200? 
Gonthier et al.l l2002t l2004li . The broad goal of these 


“Monte Carlo” (MC) simulation studies is to statistically 
reproduce the pulsar sample observed in actual surveys 
and test whether a given population model is consistent 
with the data. 

Unfortunately, the conclusions of pulsar population in¬ 
vestigations have often been conflicting. 

One long standi ng issue is that of magn e tic field de¬ 
cay. F or example. iGunn & Ostrikerl (I1970D. ILvne et al. 
(1985), iNaravan fc Ostrikei! (119901 . I Gonthier et al. 
(2002), and lGonthier et al.l ( 200411 have suggested that 
field decay occurs on short time scales ~ 2.5 — 5 Myr. 
Similar statistical studies, e.g. those of ISto llman 


1987al) , IBhattacharva et al.l (Il992f) , and ILorimer et al 


19971 ). reached the opposite conclusion that the mag¬ 


netic field of pulsars does not decay significantly dur¬ 
ing their lifetime as radio sources, implying decay time 
constants > 100 Myr. Meanwhile, theoretical argu¬ 
ments have mostly supported that field decay is unim¬ 
portant for typical neutron stars (jBavm et al.l 119691 : 
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iGoldreich fe Reiseneggerlll992i and references therein). 
Since the observed kinetic age versus characteristic age 
diagram (see, e.g., Lyne, Anderson, & Salter 1982 and 
Harrison, Lyne, & Anderson 1993, who used it to argue 
in favor of field de cay) has been shown to be an artefact 
of selection effects dLorimer et al.lH997l ). direct empirical 
evidence for field decay is also lacking. However, evolu¬ 
tion of the inclination angle between the magnetic and 
spin axes of a pulsar can also produce torque variations 
(see section 13.3.21) that are difficult to disentangle from 
evolution of the magnetic field itself on the basis of spin 
kinematics. We will thus henceforth speak more gener¬ 
ally of “torque decay” 3 . 

Ano ther matter of debate, introduced by 
IVivekanand fe Naravanl (ll981f) . concerns the “injec¬ 
tion” of a subpopulation of pulsars with birth spin 
period ~0.5 s, in addition to a population of Crab-like 
pulsars born with short period s < 1 00 ms. Evidence 

for_injection was found by IVivekanand fe Naravanl 

(il98ll) in an extension of t he “pulsa r curre nt” analysis 
proposed by IPhinnev fe Blandfordi (jl98lt ). in which 
one considers the flow of pulsars in the period-period 
derivative plane. Although the pulsar current analysis is 
model-free in the sense that it does not require explicit 
assumptions regarding the lumin o sity a nd spin-down 
laws of the pulsa rs. iLvne et all (ll985l) a rgued that 
the conclusion of IVivekanand fc Naravanl ( 1198 i f) was 
subject to considerable statistical uncertainty and 
suffered from an incomplete treatment of the selection 
effects plaguing pulsar surveys. Nevertheless, injection 
received further sup port from a number of ind epen¬ 
dent analyses (IChevalier fe E mmerinel Il986t iNaravaiJ 
Il987t iNaravan fc Ostrikeil Il990h . while simil ar studies 
maintained that the data did no t require it (Stollman 
ll987at iBhattacharva et all Il992f ). Refining the pulsar 
current analysis by restricting attention to pulsars above 
a luminosity cut-off of 10 mJy kpc 2 at 400 M Hz, for 
which significant statistics were available, iLorimer et alJ 
(Il993l) found no evidence for injection, although they 
could not exclude that it may affect the fainter end 
of the pulsar population. Injection ap pears t o have 
since lost popularity. Most recently, IVranesevic et al.l 
(l2004f ) presented a pulsar current analysis of the P arkes 
Multibeam pulsar survey dManchester et al.l 12001 ) data 
which, although it suggests that many, perhaps 40%, of 
pulsars are born with periods in the range 0.1-0.5 s, did 
not show evidence of distinct subpopulations. 

Lately, several authors have attempted to characterize 
the distribution of the birth space velocities of pulsars. 
Different analyses initially led to quan titative d isagree¬ 
ment s regarding the me a n vel o city (ILvne fc Lorimerl 
11994 lHansen fe Phinnevl Il997t ILorimer et alJ I1997D . 
Somewhat surprisingly, other workers have found 
evide nce in favor of a _ bi modal velocity distribu¬ 

tion dGordes fc ChernofflU9981 : I Arzoumanian et al.ll2002l : 


iBrisken et alJ 20031) . so that the issue of subpopulations 
of ordinary radio pulsars with different birth properties 
is not completely resolved. The origin of the two dis¬ 
tribution components has not been conclusively identi- 

3 The magnetic torque on a rotating neutron star decreases with 
time even if the magnetic field and its orientation are constant, 
owing to its decreasing rotational frequency. Here, we are referring 
to any potential decay in addition to this expected decrease. 


fied. In fact, the two components themselves remain 
to be directly exhibited in the proper-motion data and 
the different authors disagree on their modes and rela¬ 
tive weight. Analysin g a large sample of pulsar proper 
motion measurements. iHobbs et al~1 ( 2005 1 directly chal¬ 
lenged the multimodality of the distribution. 

Multiple reasons motivate us to reconsider the birth 
and evolution of isolated radio pulsars in this paper. 

First, pulsar astronomy has in the recent years seen a 
number of major advances. The Parkes ( PM) an d S win¬ 
burn e (SM) Multibeam puls ar surveys ([Manchester et aD 
1 20011: lEdwards et al.l l2001h at 1.4 GHz have approxi¬ 
mately doubled the number of known pulsars. Over 
1500 objects are now cataloged in the Australia 
Telescope Nationa l Facilit y (ATNF) Pulsa r Database 4 
([Manchester et al.1l2005l ). iGordes fc Lazid ([2002D have 
introduced a new model of the Galactic free electron 
density ( NE2001) to supersede t he previously standard 
model of iTavlor fc Gordesl (jl993l ) (TC93), providing an 
updated dispersion measure distance scale. New as¬ 
trometric measurements using interferometry provide a 
sample relatively free of brightness and distance biases 
from wh ich to estimate the space velocity distribution of 
pulsars ([Brisken et al.l[20021120031) . 

Second, in spite of poorly understood underly¬ 
ing physics, the recent simulations have become in¬ 
creasingly sophisticated. Departing from the com¬ 
mon practice of modelling the pulsar pseudo-luminosity 
(where pseudo-luminosity x distance 2 = flux density), 

I Arzoumanian et~aH (2002) have for instance introduced 
a standard-candle geometrical model of physical luminos¬ 
ity. While this luminosity model is certainly a valuable 
step toward more realistic simulations, it is somewhat 
speculative in several respects (e.g., Gaussian shape and 
angular size of the radio beams, and the relative radi¬ 
ated power contri bution of each). An i mport ant feature 
of the analysis of I Arzoumanian et alJ i (j2002t ) is the use 
of indirect information (i.e. other than proper motions) 
contained in the observed pulsar sample in the infer- 
rence of their birth velocity distribution of pulsars, re¬ 
lying on a detailed modelling of the selection effects of 
major radio surveys at 400 MHz. Since the accuracy of 
this modelling is limited by the validity of the assump¬ 
tions made regarding the intrinsic properties of the pul¬ 
sars, it is a fair question to ask whether their conclusion 
that the birth space velocity distribution of isolated pul¬ 
sars has two well defined components, corresponding to 
low (~ 90 km s _1 ) and high velocities (~ 500 km s~ 4 ) 
is really required by the data. A previous analysis by 
IGordes fc GhernoU (jl998f ) omitting a detailed treatment 
of the selection effects suggest that evidence for two com¬ 
ponents is indeed p resent in the data , but this claim has 
been disputed bv lHobbs et al.l (120051). Adoptin g a slight 
modification of the I Arzoumanian et al.l ( 20021) luminos- 
ity model to include spectral dependence. iGonthier et abl 
(12004T) claimed evidence for magnetic field decay on a 
time scale ~ 2.8 Myr. Because pulsars are ultimately 
detected through their electromagnetic fluxes, it appears 
that the conclusions of pulsar population simulations 
are highly dependent on the assumed luminosity model. 
Since there is no strong independent support for any par¬ 
ticular luminosity model, we may also consider whether 

4 http://www.atnf.csiro.au/research/pulsar/psrcat/ 
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the complexity of some of the simulations is absolutely 
needed. 

Throughout, we follow two guiding principles. First, 
because of the important uncertainties in modelling pul¬ 
sars from birth to detection, we rely on independent re¬ 
sults obtained by more direct means whenever possible. 
Second, inspired by Ockham’s razor, we strive for sim¬ 
plicity over sophistication. We fu rther attempt t o be as 
self-consistent as possible. While iGonthier et al.l (|2004 ) 
implemented the NE2001 model, for instance, they based 
their luminosity mo del and the ir birth veloc ity distribu¬ 
tion on those of lArzoumanian et ahl (12002 1, who used 
TC93. 

We begin by estimating the pulsar birth space velocity 
distribution directly from proper motion measurements 
in section [2] In section 0 we fix this distribution and 
proceed with Monte Carlo simulations of pulsar birth, 
evolution, and detection. We discuss our results and their 
implications in section [5] and conclude in section [5] 

2. THE PULSAR BIRTH SPACE VELOCITY DISTRIBUTION 
FROM PROPER MOTION MEASUREMENTS 

Pulsars have long been recognized as a high-velocity 
stellar population. Based on the correlation between the 
ages of pulsars and their distances from the Galactic 
plan e, and on the tra n sverse velocity of the Crab pul¬ 
sar, iGunn fc Ostrikerl (|l970l ) proposed that most pul¬ 
sars are born with a space velocity ~100 km s -1 . We 
now know that, in fact, some pulsars may have veloc¬ 
ity exceeding 1000 km s -1 . PSR B2224+65 in the Gui¬ 
tar Nebula, for instance, has infer red tranverse velocity 
1640 km s" 1 (|Cordes et all 119931: IChatteriee fc Cordesl 
l2002tl2004l) . Although this estimate assumes the NE2001 
dispersion-measure distance of 1.9 kpc, the observed bow 
shock is a strong indication of unusually rapid mo¬ 
tion. PSR B1508+55 has a velocity determined di¬ 
rectly from proper -motion and parallax m easurements of 
1083+go 3 km s" 1 (IChatteriee et all l2005f b These veloci¬ 
ties are to be contrasted with those ~15 km s -1 typical 
of their pur ported progenitors, massive main sequence 
stars (e.g.. lLeaueuxnl979t h 

The identified mechanisms from which neutron stars 
may acquire their high initial velocities fall in two cate¬ 
gories. The original proposal is that the velocities origi¬ 
nate from recoil in binaries th at are disrupted by a sym¬ 
metr i c supernova explosion dBlaauwl 119611 : iGott et all 
119701 : llben fc Tutukovl Il996f) . The other involves su¬ 
pernova asymmetries, which may induce “kicks” to 
nascent neutron stars. There is strong evidence that 
this second scenario plays an important role. The most 
direct evidence comes from the observed characteris¬ 
tics of binary systems containing a neutron star. In 
these systems, the spin and orbital angular momen¬ 
tum vectors are expected to be aligned prior to the 
last su pernova, o wing to tid al and mass transfer ef¬ 
fects dBliattacharva fc van den Heuvell Il99lh . A post¬ 
supernova misalignment then indicates that the super¬ 
nova must have imparted a kick to the neutron star 
out of the orbital plane, which can result only from 
an asymmetric explosion. Spin-orbit misalignment has 

in fact been detected in several bin ary system s:_in 

the PSR J0045—7 319/B-star system (lLai et al.l 119951 : 
iKaspi et al.l I1996D , in the PSR J17 40—3052/massive 
companion system (|Stairs et al.ll200,‘il) . and in the dou¬ 


ble neutron star systems PSR B1913+16 (iKrairieill 19981 : 
IWex et al.ll200l. PSR J0737-3039 and P SR B1534+12 
( Willems et al.ll2004 : iThorsett et al.ll2005l) . Direct polar¬ 
ization observations of super novae also sugg est that that 
most are asymmetric (e.g JWang et al.lf200lh . In general, 
both binary breakup and kicks due to supernova asym¬ 
metries will contribute to the observed pulsar velocities. 


As already mentio ned , several author s (e.g . 


Lvne & Lorimei 

11994 Hansen & Phinnevl 19971; 

Hartman 19971: 

iLorimer et al. 1997T Cordes & Chernofl 


1998t lArzoumanian et al.ll2002l ; IBrisken et al.ll2003l) have 


attempted to constrain the birth velocity distribution 
of pulsars. They have found mean birth velocities 
~ 300 — 500 km s _1 . All of these studies (except for 
Lorimer et al. 1997, who used the older model of 
Cordes et al. 199 1) relied on the TC93 model. As was 
demonstrated by iLvne fc Lorimei] (I1994D . adoption of a 
new free electron density model can have a considerable 
impact on derived velocities. In the absence of pulsar 
velocity studies based on the NE2001 model in the 
literature when this part of the present work was carried 
out, we reanalyze the s ample of proper motions of 
IBrisken et all ( 20021 . 12003 ). This sample consists of 34 
pulsars, 8 of which have parallax measurements. In 
general, the pulsars in the sample are fainter and more 
distant than in previous proper motion surveys, reducing 
the bias against the high-velocity objects which rapidly 
escape the solar neighborhood. Moreover, the measure¬ 
ments are typically much more precise than in most 
earlier works, which also appear to have reported overly 
optimistic error bars, most likely becaus e of a failure 
to fu lly account for ionospheric effects ((Brisken et al.l 
120031 ). For these reasons, we choose not to include other 
catalogued proper motion measurements. A caveat, 
however, is that the small sample size provides limited 
statistics. 

The basis of our formalism is identical to that of 
IBrisken et alJ (|2003f) . The Bayesian maximum likelihood 
framework has the benefit of allowing a detailed treat¬ 
ment of the uncertainties on the measured and derived 
quantities. This is important, because distance uncer¬ 
tainties, for example, are in general neither Gaussian 
nor symmetric, and ignoring this would lead to biased 
derived velocities. The consideration of non-Gaussian 
models and the calculation of credibility ranges are our 
additions. 

In section EU we describe the maximum likelihood 
formali sm u sed. The mode ls investigated are presented 
section 12.21 In section 12.31 we present and discuss the 
results. 

2.1. Method 

In this paper, we define the term birth velocity as the 
post-supernova velocity of the pulsar relative to its local 
standard of rest (LSR), including any contribution due 
to pre-supernova binary orbital motion or motion of the 
progenitor system relative to the LSR. 

Let fii = l (where l and b are the Galactic longitude 
and latitude, respectively) be the (observed) component 
of proper motion parallel to the Galactic plane, and D 
be the distance to the pulsar. Then vi = D^i cos b — A vi 
is the component of the pulsar’s transverse velocity par¬ 
allel to the plane, relative to its LSR, where A vi(D, l, b) 
is the contribution to the observed velocity due to dif- 
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ferential Galactic rotation and the motion of the Sun 
relative to its own LSR (see, e.g., iBinnev fc Merrifieldl 
119981) . We assume a flat rotation curve with circular 
velocity 225 km s -1 and a solar mot ion of 1 6.5 km s^ 1 
in the direction (l, b) — (53°, 25°) (|Mihalas fc Binnevl 
SMI). These values are consistent with those ob¬ 
tained from a more recent analysis of the kinematics 
of Galactic Cepheids fr om Hipparcos proper motions 
dFeast fc Whitelocklll997i) . The probability density func¬ 
tion (PDF) for Vi can be computed as 


pOO poo 


p(vi) = 


S{vi — [Dcosb p,i 


' 0 J — oo 


A vl }}p(pi)p(D)dp,idXl) 


The functions p(pi) and p(D) are the PDF for the com¬ 
ponent of the proper motion parallel to the plane and 
the distance to the pulsar, respectively, and reflect the 
uncertainties on their measurements. 

For pLi , the uncertainty is Gaussian, as reported by the 
observers. For D, there are two cases. If a parallax ( zu ) 
measurement is available for the pulsar, we consider the 
parallax distance and the reported Gaussian uncertainty 
on the parallax is transformed to the corresponding (non- 
Gaussian) uncertainty on the derived distance: 


p{D) 



exp 


(1/D — tu) 2 
2 ^ . 


( 2 ) 


For pulsars for which no parallax measurement is avail¬ 
able, we use the distance derived from the dispersion 
measure (DM) using the NE2001 model. To model the 
uncertainty on the distance due to imperfections of the 
free electron density model, a Gaussian uncertainty with 
variance 0.4 DM is assumed on the measured DM. Then 


p(D) = p(DM) 


d (DM) 
d D 


(3) 


The derivative d(DM)/dD is numerically approximated 
by interpolation from a table of dispersion measures as a 
function of distance computed using the NE2001 model. 
The distance step is set to 10 pc. In reality, the er¬ 
ror on the measured DM is negligible. The 40% fig¬ 
ure is approximately equal to the variance of the frac¬ 
tional errors of the model dispersion measures, with re¬ 
spect to the measured ones, for the pulsars in ATNF 
Pulsar Database with parallax measurements, assum¬ 
ing the NE2001 model and that the parallax distance 
is exact. The rationale behin d this rough model intro¬ 
duced bv iBrisken et al.l (l2002t ) is that the estimated un¬ 
certainties tend to be greater for pulsars in regions where 
the free electron density fluctuates rapidly, as expected. 
We note that out of 34 pulsars in the sample, 20 have 
Vi < 0 versus 14 with vi > 0. Under the assumption of 
isotropic velocities, an asymmetry of this magnitude has 
20% probability of occurring for our sample size and is 
not worrisome. 

We quantify the likelihood of models based on the dis¬ 
tribution of Vi for the pulsars in the sample. Although 
the sample is not restricted to young objects, dynami¬ 
cal evolution and the potential bias toward low-velocity 
objects that remain in the solar neighborhood for longer 
periods are not expected to significantly affect our anal¬ 
ysis. In fact, we verified that the velocity components 
parallel to the Galactic plane of the old pulsars are, on 
average, close to their birth values. Specifically, we used 


the computer model described in section [3] to compare 
the averages of the Cartesian velocity components par¬ 
allel to the Galactic plane after evolution of the pulsars 
to the averages of their birth values. The differences 
were ~ 5%, which is less than the uncertainties on the 
model parameters derived here (c.f. section m - Be¬ 
cause younger pulsars appear generally easier to detect 
(see section l4~4l) . the actual effect in the observed sam¬ 
ple is presumably even smaller. The bias due to high- 
velocity pulsars rapidly moving away from the Galactic 
plane (and hence from the sampled volume) is mitigated 
because pulsars with a large velocity component perpen¬ 
dicular to the plane do not necessarily have correspond¬ 
ingly large components parallel to it. It is implicitly as¬ 
sumed that the birth velocities of the pulsars are inde¬ 
pendent of their other characteristics. 

The likelihood for a model A4, depending on the (pos¬ 
sibly vectorial) parameter 9, is 

/ OO 

Ppsr(vi)p M (9)(vi)dvi. (4) 

-o° 

We maximize this function with respect to 9 to find the 
maximum likelihood value of the parameter. The evi¬ 
dence for the model is 

£(M) = J £(M(9)) Po (9)d9, (5) 

where po(9) is the prior PDF for the parameter and 
the integration is over the entire parameter space, i.e. 
the set of values that 9 may take. The evidences al¬ 
low comparison of different models. The “odds ratios” 
£j /£j may be thoug ht of as rel ative probabilities of mod¬ 
els i and j (ID ’ Agostinil 120031) . In all cases, flat priors 
(po = constant) are assumed over a reasonable param¬ 
eter space. The exact evidence for each model depends 
on this (somewhat arbitrary) choice of priors, though this 
dependence vanishes in the limit of a large data set. We 
quote (normalized) evidences to two digits to illustrate 
the differences between models with similar evidences, 
although only the first digit is really significant and the 
numerical values should not be overinterpreted. 

For each model, we provide an uncertainty estimate 
for the maximum likelihood parameter. Let 9ml be 
the parameter of maximum likelihood. We solve for 
9min ^ 9ml ^ 9 rnax so that 

Iei:^M0))d9 C(M(9))d9 

5T 00 C(M(9))d9 fZ£(M(0))d9 ’ 1 j 

where C is the “probability content” ([Gregory fc Loredol 
Il992t) . We consider lcr errors, for which C=34.15%. If 9 
is a vector, as in the case of the Gaussian model below, we 
treat each component as above, fixing the others to the 
maximum likelihood values. The intervals [9 m i„, 9 max \ 
are known as “credibility ranges”. 

2.2. Models 

We investigate several functional forms for the pulsar 
birth velocity distribution. We describe them in this sec¬ 
tion. 


2.2.1. Gaussian Model 
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The first model that we consider is a two-component 
Gaussian: 

. . wi ( vf \ 1 - wi ( vf \ 

pivi) = 75sr exp + v^ exp 

(7) 

where w\ is the fraction of pulsars in the low-velocity 
component, and (T\ and a 2 are the 1-D velocity disper¬ 
sions for the low- and high-velocity components, respec¬ 
tively. Gaussian models are natural to investigate, as 
they frequently arise as a consequence of the central 
limit theorem. The possib ility of two components i s 
motivated by the s t udies of iCordes & Chernoll ( H)9|f) , 
lArzoumanian et ah! (|2002f) . and iBrisken et all ( 200311 . 
who all considered such models. In the case w\ = 1, this 
model reduces to an ordinary Gaussian. The correspond¬ 
ing three-dimensional velocity distribution is Maxwellian 
and has mean (v^ d) — i/8/7r[wicri + (1 — 



Transverse Velocity Projected onto the Galactic Plane (km s 

Fig. 1.— Histogram of the tranverse velocity components 
parallel to the Galactic plane for the pulsars used in our anal¬ 
ysis of the birth velocity distribution. For this histogram, the 
velocity is calculated directly using the most probable distance 
(from parallax if available; otherwise derived from the DM us¬ 
ing the NE2001 model) and proper motion of each pulsar, with¬ 
out taking uncertainties into account. The overlaid curves illus¬ 
trate the maximum likelihood probability density function, fit¬ 
ted to the data with uncertainties as described in section 12.11 
for each of the investigated models: single-component Gaussian 
(dashed), two-component Gaussian (dashed-do tted), exponentia l 
(solid), Lorentzian (dashed-triple-dotted), and IPaczvnskH (119900 
(dotted). 


As illustrated by Figure [TJ while 27 of the 34 
pulsars in the sample have p;| < 300 km s -1 , the 
distribution extends beyond |u;| = 1000 km s -1 , with 
pz| = 1340 km s -1 for PSR B2011+38. Here, the ve¬ 
locity components are calculated using the most proba¬ 
ble distance and proper motion for each pulsar, without 
taking the uncertainties into account. Since Gaussian 
functions have rapidly decaying tails, we do not expect 
a single-component Gaussian to describe the data well. 
Rather, we expect a two-component Gaussian with ex¬ 
tended tails to be favored in order to accommodate the 
highest-velocity objects. Such a two-component Gaus¬ 
sian model requires three free parameters. It also re¬ 
sults in a bimodal three-di mensional velocity distribu¬ 
tion (see, e.g., Figure 3 of I Arzoumanian et all 120021 . a 
feature which, if real, requires an astrophysical expla¬ 
nation. To investigate whether the complexity of the 
two-component Gaussian model is justified, we consider 


a number of alternative single-parameter models, with 
heavier tails than Gaussian functions. 


2.2.2. Alternative Single-Parameter Models 

The (double-sided) exponential model has functional 
form 

»W = 2;yexp(-M), (8) 

where ( vi ) is the mean absolute value of the velocity 
component. The corresponding three-dimensional dis¬ 
tribution is difficult to derive analytically, but its mean 
is easily estimated numerically by simulation of velocity 
vectors. The Lorentzian model is defined by 

pM = Mvf+HWHM‘Y (9) 


where HWHM is the half width at half maximum. The 
first moment of the corresponding three-dimensional dis¬ 
tribution in this case does not exist and the mean is thus 
un defined . Finally, we consider a model introduced by 
iPaczvnsO (ll990f) : 


P(vi) 


[1 + (' vi/v*)} In [(1 + (vi/v*) 2 )/(vi/v*) 2 ] - 1 
7r[l + (vi/v*) 2 ] 


( 10 ) 

This vi distribution corresponds to the more mathemat¬ 
ically appealing three-dimensional distribution 


p(v 3 n) 


4 

7TU*[1 + ( V 3D /V *) 2 ] 2 


( 11 ) 


The parameter u* is related to the three-dimensional 
mean by (i> 3 d) = 2u*/7t. 


2.3. Results and Discussion 

The results of the maximum likelihood optimization 
are summarized in Table [D We also list the derived 
mean three-dimensional velocity and the fraction of pul¬ 
sars with birth velocity < 60 km s -1 , the typical cen¬ 
tral escape velocity for the most m assive G alactic glob¬ 
ular clusters, such as Tucanae 47 (IWebbinkl Il985h . for 
each model. As a check on our calculations, we also 
give the Kolmogorov-Smirnov (K-S) P-values Pr-S ( see > 
e.g., Press et al. 1992 and section EU) to compare the 
observed distribution of vi and the maximum likelihood 
distribution for each model. For the K-S P-values, only 
the most probable values of the proper motion and dis¬ 
tance are used to calculate the velocity component vi, 
i.e. the uncertainties are ignored, and we only expect 
the results to roughly validate those of the more detailed 
maximum likelihood analysis. A further benefit of the 
K-S P-values is that they provide an absolute measure 
of the goodness of fit, whereas the maximum likelihood 
analysis only provides relative figures. 

With an evidence ~ 10~ 3 that of the other models, 
and a definitely lower K-S P-value, the single-component 
Gaussian model is clearly inadequate. The four other 
models all have identical evidence, up to the significant 
digit, and we cannot distinguish among them. None 
of them can be rejected at a significant level of con¬ 
fidence based on the K-S test either. In particular, 
we cannot justify the complexity associated with the 
two-component Gaussian model. Nonetheless, it is 
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TABLE 1 

Maximum Likelihood (ML) Birth Velocity Models 


Model 

ML Parameters 
(km s -1 ) c 

£ a 

Searched Range 
(km s -1 ) 

(l>3 d) 
(km s _1 ) 

P K-S 

/eo b 

% 

1-Gaussian 

a = 

290i“ 

9.7 x 10“ 4 

[100, 500] 

460l 5 5 ° 

0.045 

0.2 

2-Gaussian 

W\ = 

0-90iao° 

1.0 

[0.05, 1.00] 

350t?° 0 

0.442 

1.2 


(Ti = 

= 


[50, 300] 





= 



[300, 1000] 

ssolgg 



Exponential 

v = 

180t” 

1.0 

[50, 300] 

0.377 

1.4 

Lorentzian 

HWHM = 100I,S 

1.3 

[50, 350] 

d 

0.245 

2.0 

Paczynski 

V* = 

560t“° 

1.5 

[200, 900] 

360 

0.153 

13.5 


a The quoted evidences have been normalized so that £ = 1 for the exponential model. b Fraction 
of pulsars with birth velocity < 60 km s _1 , the ty pical central esc ape velocity of the most massive 
Galactic globular clusters, such as 47 Tucanae dWebbinki 1198 S). It is 2.8% for the model of 
lArzoumanian et alj 1120021 ). c Except for the dimensionless fraction w\.^ The expectation value of 
v 3 u is undefined. 


interesting to note that the maximum likelihood pa¬ 
rameters for the two-component Gaussian model place 
u;i=0.901o'J 3 of pulsars in the low-velocity compo¬ 
nent, with 1-D dispersion <ti= 160 P 3 q km s _1 , and the 
rest in the high-velocity component, with dispersion 
CT 2 =780ti4o km s 1 (implying a mean three-dimensional 
velocity ( 113 . 0 } =350 P i4 O km s -1 ). ICordes fe ChernoU 
(119981) arrived at a similar result, with best-fit two- 
component Gaussian model parameters w\ = 0 . 86 ( l lo'i 2 ) 
0-1=175^24 km s _1 , and cr 2 =700(ti32 km s _1 (with cor¬ 
responding {v 3D )=395±Iq 5 km 1 ). lArzoumanian et ali 
( 120021 ) found maximum likelihood for u>i= 0 . 4 P q' 2 , 
<Ti=90(ti5 km s- 1 , and cr 2 =500(ti5o km s _1 (with 
corresponding (h3d)=540 P 24o km s _1 ). Fi¬ 

nally, iBrisken et al.l (120031) obtained wii=0.20, 
<ti=99 km s _1 , and 02=294 km s _1 (with corresponding 
(^3 d)= 407 km s _1 ). Here, the quoted uncertainties 
on the model parameters w±, or, and are at the 
lo- level in each case. We estimated the uncertainties 
on the derived (1130) by calculating the minimum and 
maximum values allowed when each of the parameters 
is confined in its la region. While our derived mean 
three-dimensional velocity is consistent with each of 
these results within lo - , the uncertainties are large due 
to the multiple model parameters. Moreover, the fact 
that the values of the three model parameters describing 
the two components of the velocity distribution vary 
wildly (especially the fraction of low-velocity pulsars) 
suggests that the true velocity distribution does not 
have two well defined components. This supports the 
hypothesis that previous authors were pushed toward a 
more complex two-component Gaussian model because 
of the non-Gaussianity of the true distribution, for 
instance to accommodate the more extreme objects, 
rather than because of a genuine bi-Gaussianity. 

Except for the single-component Gaussian model, for 
which the derived mean three-dimensional velocity is 
overestimated because of the highest-velocity pulsars, 
the maximum likelihood mean three-dimensional veloc¬ 
ity for each of the models we have investigated, which 
varies from 350 km s -1 to 380 km s -1 , is smaller than 
most previous estimates (e.g., 450 km s -1 for Lyne & 
Lorimer 1994, ~500 km s ” 1 for Lorimer, Bailes, & Har¬ 


rison 1997, and the results discussed in the previous 
paragraph). This is consistent with the distances esti¬ 
mated using the NE2001 model being in general smaller 
than those obtained using TC93, which is the case 
for all the pulsars with distances derived from dis¬ 
persion measures retained in our analysis, except for 
PSRs B1534+12 and B1541+09. A notable exception is 
lHansen &; Phinnevl (|l997l) . who found a mean birth ve¬ 
locity ~ 250 — 300 km s^ 1 . The origin of the discrepancy 
is not clear, but is pos sibly attributable to pulsar weight 
adjustments made by lHansen fe Phinnevl (I1997T) to at¬ 
tempt to account for Malmquist-like biases (B. Hansen 
2005, private communication). 

There are two main diff erences betwee n our con¬ 
clusions and those of IBrisken et al.l ( 2003 1, who used 
the same pulsar sample and performed the same ba¬ 
sic analysis. First, our best-fit velocity models, exclud¬ 
ing the single-component Gaussian one, have lower (al¬ 
beit consistent within ler) mean three-dimensional ve¬ 
locities (~ 350 — 380 km s -1 compared to 407 km s” 1 ). 
This is expected, as just discussed, since we have 
used the NE2001 mode l for the G alactic free electron 
dens ity, while IBrisken et, aD (1200.31 ) used TC93. Sec¬ 
ond, IBrisken et al.l (20031 ) favor a two-component Gaus¬ 
sian functional form, while we prefer alternative single¬ 
parameter model s. We understand this to be simply 
a result of IBrisken et al.l ( 2003 ) not investigating non- 
Gaussian single-parameters models, for had we only con¬ 
sidered Gaussian models as they did, we too would have 
been led to favor a distribution with two components (see 
Tabic [T]). 

After we had completed our analysis, iHobbs et all 
(200.5) published a pulsar birth velocity distribution 
inferred from a sample of 73 pulsars with measured 
proper motions (most of which from pulse timing) and 
characteristic ages r c = P/2P <3 Myr using NE2001. 
They did not, however, treat measurement uncertain¬ 
ties. They found that the data are well described 
by a three-dimensional Maxwellian distribution with 
one-dimensional standard deviation a = 265 km s — . 
This standard deviation is consistent with the one 
we obtained for our single-parameter Gaussian model 
(a = 290±|g km s *), but the result is at odds with the 
fact that w e disfav o r a single-p arameter Gaussian func¬ 
tional form. IHobbs et al.l (120051) did not investigate alter- 
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native single-parameter models and we note that PSRs 
B2011+38 and B2224+65 (both with inferred transverse 
velocity > 1600 km s -1 ), which they have included in 
their analysis, have vanishingly small probability of oc¬ 
curring under the hypothesis that the velocities are dis¬ 
tributed according to their best-fit Maxwellian distribu¬ 
tion 5 . Rather than directly fitting a model to the ob¬ 
served one- or two-dimensional distribut ions of ve l ocities 
projected on the plane of the sky, iHobbs et al.l (120051 ) 
used a deconvolution algorithm to estimate the corre¬ 
sponding three-dimensional velocity distribution, under 
the assumption of isotropy of the velocity vectors. They 
then fitted a model to the deconvolved data. The fact 
that PSRs B2011+38 and B2224+65 are inconsistent 
with their best-fit model, despite the fit to the decon¬ 
volved data being statistically good (reduced y 2 = 0.6) 
suggests that some information was lost in the deconvo¬ 
lution process and that the technique may not be sensi¬ 
tive to subtle tail behavior. It is true, though, that the 
authors tested their deconvolution algorithm and veri¬ 
fied that it reproduced mock three-dimensional distribu¬ 
tions well, at least to a degree sufficient to distinguish 
between unimodal and bimodal distributions. We do 
no t necessar ily c onside r our results discrepant with those 
of IHobbs et al.1 (120051 ). as the difference appears to lie 
mainly in whether a few extreme objects in the sample 
are given weight. The resolution of this dilemma will 
require reliable confirmation or rebuttal of the measure¬ 
ments of the high velocities and will depend on whether 
such obj ects con tinue to be discovered. 

lHansen fe Phinnevl (119971) have remarked that there 
is significant uncertainty regarding the form of the pul¬ 
sar birth velocity distribution. They demonstrated that, 
due to projection effects, a three-dimensional kick ve¬ 
locity distribution consisting of two delta functions (lo¬ 
cated at 250 km s _1 and 1000 km s , with weight 0.8 
and 0.2, respectively), was equally consistent with the 
observed proper motions as their preferred, qualitatively 
very different, Maxwe llian with one- dimensional disper¬ 
sion a =190 km s -1 . iLorimer et al.l ( 1997 ) also showed 
that a variety of birth velocity distributions were consis¬ 
tent with the data. We continue to find that the exact 
shape of the birth velocity distribution is poorly con¬ 
strained, as illustrated by our inability to quantitatively 
discriminate between four different models. Reliable de¬ 
termination of the correct shape of the kick velocity from 
proper motion measurements will require study of a con¬ 
siderably larger sample of pulsars. Perhaps our non¬ 
detection of two well defined components is due to the 
size of the analyzed sample, but given the lack of com¬ 
pelling evidence for a multimodal distribution, we follow 
Ockham’s razor in favoring unimodal, single-parameter 
models. 

The cumulative three-dimensional velocity distribution 

5 Consider n three-dimensional velocities independently drawn 
from a Maxwellian distribution with standard deviation a for the 
corresponding one-dimensional Gaussian distribution: = 

exp (—Ug£j/2(7 2 ) (v^d >0). Each velocity has prob¬ 
ability p = p(v 3 o)dv 3 £> of exceeding the value v. For in¬ 

dependently drawn velocities, the total number E of veloci¬ 
ties exceeding v is binomially distributed with parameters n 
and p. The probability that E is at least 2 is thus given by 
1 — (1 — p) n_1 (l + (n — l)p). For a = 265 km s -1 , n = 73, and 
v = 1600 km s -1 , we obtain the value ~ 10“ n . 



Fig. 2.— Cumulative distribution function for the three- 
dimensional birth velocity for each of the maximum likeli¬ 
hood models investigated: single-component Gaussian (dashed), 
two-component Gaussian (dashed-do tted), exp onential (solid), 
Lorentzian (dashed-triple-d otted), and Paczynski (199(f) (dotted). 
The distribution inferred b v I A rzo: ] m a n i a rTTiaTTl 200T ) is indicated 
by the thick solid curve. 

for the maximum likelihood parameters for each model 
considered here is shown in Figure [2j The distributions 
are consistent with the observed retention rates of neu- 
tron stars i n globular clusters. In fact, the results of 
I Ivanova et al.l (120051) on the dynamical formation and 
evolution of binaries containing neutron stars in dense 
clusters suggest that concerns regarding the fraction of 
neutron stars retained in globular clusters are not justi¬ 
fied (the so-calle d “retentio n problem”; for a discussion, 
see for example lPfahl et al.ll2002f ). This is thanks to neu¬ 
tron stars in clusters being very effectively recycled and 
the total number of neutron stars in clusters being much 
smaller than previously thought (F. Rasio 2005, private 
communication). In particular, in light of a deep Chan¬ 
dra X-Ray Observatory observation, the globular cluster 
47 Tuca nae is n ow estimated to contain ~25 millisecond 
pulsars dHeinke et al.ll2005l ), compared to a pr evious esti¬ 
mate b ased on radio observations of >200 ([Camilo et al.1 
l2000bD . 

In most of the remainder of this paper, we adopt the 
single-parameter model with the smallest fractional un¬ 
certainty (15%): the exponential one. The corresponding 
three-dimensional mean velocity is 380lgg km s -1 . We 
consider how our subsequent results would be affected 
if we had chosen one of the other models with similar 
evidence in section roi 

3. POPULATION SYNTHESIS 

We now proceed with a more general pulsar population 
synthesis. We consider only isolated, non-recycled radio 
pulsars in the Milky Way field. 

3.1. Galactic Model 

As the first step, we must define a model for the 
galaxy in which the synthetic pulsars will be placed. 
We model the spiral arm structure of the Milky Way, 
its large-scale gravitational potential, and its distribu¬ 
tion of free electrons. We introduce a Galactocentric 
system of coordinates. The origin is defined to be the 
Galactic Center (GC). The x, y, and 0 axes, respectively 
parallel to ( l,b ) = (90°, 0°), (180°, 0°), (0°, 90°), form a 
right-handed Cartesian frame. As usual, r = - \Jx 2 + y 2 , 
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TABLE 2 

Spiral Arm Parameters 


Arm Number 

Name 

k 

(rad) 

r 0 

(kpc) 

e 0 

(rad) 

i 

Norma 

4.25 

3.48 

1.57 

2 

Carina-Sagittarius 

4.25 

3.48 

4.71 

3 

Perseus 

4.89 

4.90 

4.09 

4 

Crux-Scutum 

4.89 

4.90 

0.95 


Note. — The nume rical values differ from those given by 
IWainscoat et afl 1119921 ! due to differences between the coor¬ 
dinate systems used in that paper and in this one. 


6 = arctan (y/x), and x = (x, y , z). 

3.1.1. Spiral Arm Structure 

Isolated neutron stars are thought to be formed in core¬ 
collapse supernova events. These mark the deaths of 
massive, short-lived Population I stars associated with 
the arms of spiral galaxies. It is thus expected that the 
birth sites of isolated pulsars will be hig hly co r related 
with Galactic spiral arms. In fact. iKramer et ahl (|2003al 
see especially their Figure 5) note that following the PM 
survey, the Galactic spiral structure is now clearly vis¬ 
ible in the distribution in Galactic longitudes of young 
pulsars, making it necessary for any realistic pulsar pop¬ 
ulation synthesis study to model this spiral structure. 
This has not previously been done in any of the major 
published works. 

We model the spiral structure of the spatial distribu¬ 
tion of the pulsar progenitors by four major arm cen¬ 
troids whose loci are described analytically by equations 
of the form 

9{r) = k In(r/r 0 ) + 6 0 . (12) 

The values for the parameters for each arm, given in 
Table [2] are taken from Wainscoat et al. (1992) and 
are consistent with those of the NE2001 model. We do 
not, however, model the “local arm” nor the perturba¬ 
tions to the spirals included NE2001. We adopt a Sun- 
GC distance of Rq = 8.5 kpc. This distance is consis- 
tent with NE2001, Hipparc os proper motions of Cepheids 
(jFeast fc Whitelocki ll99 7j. and a single-step measure- 
rnent using red clump stars dPaczvnski fe Staneklll998ll . 


3.1.2. Gravitational Potential 

We adopt the m odific ation by Kuiiken fe Gilmorel 
ljl989f) of the lCarlberg fe Innanenl (119871 ) fit of the Galac¬ 
tic gravitational potential. This model consists of a disc- 
halo component, a bulge component, and a nucleus com¬ 
ponent: 

<t>G{r, z) = (j)dh(r , z) + 4>b(r ) + 0„(r), (13) 


where 
<t>dh{r,z ) 


and 


-GM dh 

+ E-=i A v^ + 7^) 2 + b 2 dh + r* 

(14) 


4>b,n(r) 


-GM bi n 

\J b ln + r2 


(15) 


The numerical values of the constants are given in Table 

m 



x (kpc) 

Fig. 3. — Example of a simulated initial distribution of pulsars 
in the xy plane. Each point represents the birth position of a 
pulsar projected onto the Galactic plane. The location of the Sun 
is indicated by the symbol 0. Solid lines trace the spiral arm 
centroids. Dashed lines are tangent to arms in the vicinity of the 
Sun. 


TABLE 3 

Parameters for the Galactic Gravitational Potential 
Model 


Constant 

Disc-Halo (dh) 

Bulge (b) 

Nucleus (n) 

M 

1.45 x 10 11 M 0 

9.3 x 10 9 M 0 

1.0 x 10 10 M 0 

di 

0.4 




0.5 



P 3 

0.1 



hi 

0.325 kpc 



h2 

0.090 kpc 



h3 

0.125 kpc 



a G 

2.4 kpc 



b 

5.5 kpc 

0.25 kpc 

1.5 kpc 


Note. — Values taken from IKuiiken fe Gilmorel IU989T ). 


3.1.3. Free Electron Density 

The propagation of radio pulses from the pulsars to us 
is affected by the intervening interstellar medium. First, 
it introduces a frequency dependence of the group veloc¬ 
ity of the radio waves which causes the observed pulses 
to be dispersed by an amount proportional to the inte¬ 
grated column density of free electrons (a quantity known 
as the dispersion measure, DM). Second, the density fluc¬ 
tuations of the free electrons scatter the radio waves and 
the resulting multipath propagation broadens the pulses. 
These two effects hamper the detectability of the pulsars 
and must be modelled. To do so, we use the software 
implementation of the NE2001 mode l of Galactic free 
electron density (jGordes fc Lazioll2002l ). which, given the 
position of a synthetic pulsar in the Galaxy and the ob¬ 
serving frequency, computes its modelled DM and pulse 
broadening scattering time (r scaM ). To reduce the com¬ 
putational cost, we have changed the default integration 
distance step of 10 pc to 50 pc. 
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3.2. Pulsar Birth Properties 

The modelled Galaxy is populated by synthetic pul¬ 
sars. We describe how birth properties are assigned to 
them. 

3.2.1. Spatial and Kinematic Distributions 


The birth spin period Po of each pulsar is chosen from a 
normal distribution with mean (Po) and standard devia¬ 
tion CTp 0 . Negative spin periods are rejected and redrawn. 
The equatorial surface magnetic field B is chosen from 
a normal distribution with mean (log B) and standard 
deviation ai og B in the logarithm to the base 10. 


The birth locations of the Monte Carlo pulsars are 
specified by a radial (r) distribution and a vertical (z) 
distribution. 

The spiral structure is realized by choosing the birth 
locations so that their projections lie on arms (assuming 
an equal pulsar birthrate in each arm) and subsequently 
altering them to simulate a spread about the arm cen¬ 
troids. More precisely, an arm number, 1 < i < 4, is first 
randomly 6 chosen. A distance r raw from the GC is then 
chosen according to the specified radial distribution (see 
below) and the corresponding polar angle 9 raw is cal¬ 
culated so that the pulsar lies on the centroid of arm 
i. To avoid artificial features near the GC, the distri¬ 
bution is blurred by applying a correction of magnitude 
9 CO rr exp(—0.35r raw /kpc), where 9 corr is randomly cho¬ 
sen in the interval [0, 2n) rad, to the polar angle of each 
pulsar. Finally, to spread the pulsars about the spiral 
centroids, the (x, y) coordinates of each pulsar on the 
Galactic plane are further altered by translating it by 
a distance r corr drawn from a normal distribution cen¬ 
tered at zero with standard deviation 0.07r rau ,, without 
preference with respect to direction. These corrections 
are somewhat arbitrary, chosen so that they produce a 
reasonably natural distribution. Figure [3] illustrates a 
resulting simulated birth distribution in the xy plane. 

For the radial distribution of the puls ar pro geni tors, we 
adopt the functional form suggested bv lYusifov fe Kiicukl 
(120041 ) following an analysis based on the NE2001 model 
and including data from the PM and SM surveys: 


p(r) = A 


( r + Ri 

\ R® + Ri 


a 

exp 



r-Rp V 
Rq + Ri) 


(16) 


where p(r) is the surface density at Galactocentric ra¬ 
dius r; R \, a and b are model parameter s ; and A is a 
normalization constant. lYusifov fc Kilciik (j2004t ) found 
the best-fit values a = 1.64 ± 0.11, b = 4.01 ± 0.24, 
and R\ = 0.55 ± 0.10 kpc. Within the limits of the esti¬ 
mated errors, this radi al d istribu tion is nearly the same 
as the on e derived by iLorim e il (120031) using a d i fferent 
method. lYusifov fe Kiiciikl (|2004D and iLorimeil (|2003l ) 
actually investigated the distribution of evolved pulsars, 
rather than that of their progenitors, and it is a priori 
unclear whether it is justified to adopt this model for the 
birth distribution. We address this question a posteriori 
in section 14.1.11 

The distance from the Galactic plane is simply chosen 
from a exponential distribution with specified mean, (zo), 
and is assigned a random sign. 

To model the birth space velocities of the pulsars, we 
adopt the exponential velocity distribution determined in 
section[2j We examine the importance of this assumption 
in section ITTfl 


3.2.2. Spin Period and Magnetic Field 


6 In this paper, the unqualified term “random” implies an uni¬ 
form distribution. 


3.2.3. Radio Luminosity 

Despite nearly four decades of study, the radio emission 
mechanism and the geometry of pulsar emission beams 
remain poorly understood. Nonetheless, we must model 
the pulsar luminosity in our simulations in order to de¬ 
termine which objects are detectable. 

To avoid modelling potential apparent luminosity vari¬ 
ations due to complex emission beam and viewing ge¬ 
ometries, each pulsar is assigned a (pseudo-)hmhnosity 
L defined such that L = SD 2 , where S is the pulsar’s 
flux density at the Earth, instead of a true physical lu¬ 
minosity. 

Given the lack of well defined correlation of pul¬ 
sar radio luminosities with other properties in the ob¬ 
served pulsar sample, a model that must be tested is 
one in which the luminosities are assumed indepen¬ 
dent of all other properties. This mo del will be re¬ 
ferred to as the “random” model. iLvne et ~aD (119981 ) 
investigated the luminosity distribution of pulsars at 
400 MHz. For L 400 > 10 mJy kpc 2 , they found that 
-^(A4c>o)dT4oo ^400^400- As L 400 —> 0 , this law must 
break down, for otherwise the integral giving the total 
number of pulsars in the Galaxy is divergent. We model 
a flattening of the luminosity distribution by a second 
power-law component for luminosities below a certain 
turn-over point L to , with continuity enforced at L to , and 
a low-luminosity cut-off Li ow ensuring convergence: 

( L ai for L e [L low ,L to ) 
p(L) cx < L° 2 for L S [L to ,oo) . (17) 

[ 0 otherwise 

Assuming that the spectral index is independent of the 
luminosity at a particular frequency, the distribution of 
luminosities at 1.4 GHz should have the same shape as 
at 400 MHz, except for bein g shifted toward lower lu¬ 
minosities. From Figure 8 of ILvne et al.l (11998D. we es¬ 
timat e that for spectral indexes ~ —1.8 ( Maron et al.l 
12000D . L to = 2 mJy kpc 2 , ai = —19/15, and aq = —2 
at 1.4 GHz. As discussed in section 14.4.11 our conclu¬ 
sions are unaffected by this crude approximation. We 
set Li ow = 0.1 mJy kpc 2 , roughly corresponding to the 
faintest observed pulsars. 

Since the radiated energy is thought to be derived 
from the loss of rotational kinetic energy due to mag¬ 
netic braking, another natural model to test is one 
in which the radio luminosity is related t o pul sar’s 
period and period derivative. ILvne et all ( 1975 ) ar¬ 
gued that the pulsar luminosity has a power-law depen- 
de nce on P and P. IVivekanand fc Naravanl (11981 ) and 
iProszvnski fe Przvbvcienl (119841) reached similar conclu¬ 
sions, and this has often been assumed in later works. 
The correct power-law exponents, though, are still un¬ 
determined. Following the common practice, we also in¬ 
vestigate a model in which L is proportial to a power 
law in P and P. We dither the standard candle luminos¬ 
ity by adding to it a correction in the logarithm. This 
presumably accounts for both physical variations about 
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TABLE 4 

Measured Pulsar Braking Indexes 


Pulsar Name 

Braking Index 

References 

B0531+21 (Crab) 

2.51±0.01 

1 , 2 

B0540-69 

2.140±0.009 

3 , 4,5 

B0833-45 (Vela) 

1.4±0.2 

6 

J1119—6127 

2.91T0.05 

7 

B1509-58 

2.839±0.003 

8, 9 


References. ■ Cl I ILvne et all !ll 988fl: (2 ) 

ILvne et all 1119931: 131 ILivingstone"et all 11 2005a ); 
f4)IGouifes^t*ET[j [lj92l L f 5jjDeete^e^5Ti^i9£ ); 

(6) U^6||; (7) [Camilo^etjiyJ^^^ b 

f ^^ JRaspr’etaL ll 1 9941) ; f9) ILivmgstone"et**aLl 
12005bll 

the modelled luminosity and observed variations due to 
differing viewing geometries: 

logL = log (^L 0 P £p P 15 ep ^J + L corr , (18) 



where P is in s, Pi 5 is in 10 -15 s s -1 , and L corr is cho¬ 
sen from a zero-centered normal distribution with stan¬ 
dard deviation <JL corr - The required dither about the 
luminosity law is also likely to be partially due to er¬ 
rors in the distance measurements to actual pulsars (see 
section [4.1.ID . which we do not treat explicitly here. In 
this “P — P power law” model, luminosities are time- 
dependent. 

Since we simulate only radio surveys at 1.4 GHz (see 
section 13.51) . we consider only luminosities at that fre¬ 
quency. 


3.3. Pulsar Evolution 

After birth properties are chosen, the synthetic pulsars 
are evolved in time. 

3.3.1. Spatial Evolution 

The empty space between stars dominates the vol¬ 
ume of the Galactic field , so that enc ounters are ex¬ 
tremely rare dBinnev fe Tremainel I1987T) . As a conse¬ 
quence, the orbits of pulsars are determined by the large- 
scale, smooth Galactic gravitational potential according 
to the equation 

x — - V0 G , (19) 

where (j>G is given by equation [TB] The orbits are 
computed numerically using the adaptive time-stepping 
LSODA solver of the Lawrence Livermore National Labo¬ 
ratory ODEPACK 7 collection of ordinary differential equa¬ 
tion solvers. 


3.3.2. Rotational Evolution 


The losses of rotational kinetic energy due to radiation 
by pulsars are inevitably associated with losses of angu¬ 
lar momentum and the lengthening of the spin period. 
These losses are generally assumed to be dominated by 
magnetic dipole braking. In the case of a perfect dipole 
rotating in a vacuum, 


PP = 


/8t r 2 P 6 \ 


B 2 sin 2 x, 


( 20 ) 


where R is the radius of the star, I is its moment of iner¬ 
tia, c is the speed of light, and x is the angle determined 


7 http://www.llnl.gov/CASC/odepack/ 


Period (s) 

Fig. 4. — P — P diagram for the real pulsars retained in our 
analysis. Each point indicates the period and period derivative 
of a pulsar. Contours of constant surface magnetic field and of 
constant characteristic age are indicated by dashed and dotted- 
dashed lines, respectively. The thick solid line marks the modelled 
death line. 


by the magnetic and spin axes (|Ostriker fe Gunnl fl9691. 
Despite the extreme gravitational binding energies of 
electrons and protons at the surface of a neutron star, 
the magnetic field of a pulsar is so strong (with a typical 
value ~ 10 12 G) that a surface charge layer cannot be in 
dynamical equilibrium, causing the magnetosphere to be 
filled with plasma. A torque such that PP oc R 6 B 2 /I is 
then exerted on the rotati ng pulsar even if the spin and 
magnetic axes are aligned ([Goldreich fe Julianlll969t h so 
that equation [20] with the factor sin 2 x is invalid in the 
limit x —* 0. This suggests that even orthogonal rota¬ 
tors are fairly well approximated by setting sin 2 \ ~ 1 in 
equation 1201 Realistic pulsar spin down in a plasma is 
still poorly unde rsto od and is being actively researched 
(see, e.g.. ISpitkovskvl[2004f ). For simplicity, we assume 
that the pulsar spin down is determined by equation 
l20l with sin 2 x = 1 and we adopt the canonical values 
R = 10 6 cm and / = 10 45 g cm 2 (jTavlor fc Manchester ! 
H97l . Any torque scatter due to the distribution of the 
magnetic inclination angles, if in fact significant, can be 
thought of being absorbed in the distribution of the B 2 
factor, although we will continue to refer to B simply as 
the magnetic field to alleviate the text. 

These assumptions have the benefit that the magnetic 
field B of a synthetic pulsar is equal at all times to the 
magnetic field of a real pulsar with the same P and P in¬ 
ferred from pulse timing using the conventional formula 

Burning = 3.2 x 10 19 PP/s G. They also imply a brak¬ 
ing index n = 3, where P oc P 2 ~ n . The braking indexes 
of several pulsars have been measured and are in fact 
in the range 1.4—2.91 (see Table SI. We experimented 
with normal distributions of braking indexes in our sim¬ 
ulations. The results were weakly affected and it was 
not possible to constrain the correct mean and standard 
deviation of the underlying distribution. Again for sim- 
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plicity, we chose to keep n = 3 for all synthetic pulsars. 
We discuss the impact of this choice on our results in 
section 14.61 

Integrating equation [20] shows that, assuming a con¬ 
stant magnetic field, the spin period of a pulsar of age t 
given by 

p = V /p » + < 21 > 

In the case of the P — P power-law model, the evolved 
luminosity is calculated according to the evolved period 
and period derivative. 

3.3.3. Radio Emission Shut-Off 

The production of electron-positron pairs, as a source 
of particles to accelerate, has been assumed essential for 
radio pulsar emission. The polar caps, located just above 
the star’s surface at the magnetic poles, are the presumed 
sites from which radio waves originate. In most polar- 
cap models, radio emission ceases when the potential 
drop required for pair production exceeds the maximum 
which can be achi eved in the r otating pulsa r’ s magne¬ 
tosphere (ISturrocklll97ll: [Ruderman fc Sutherlandlll975l : 
I Chen fe Rudermanl 1199311 . The precise time at which 

the turn-off occurs depends on both the structure and 
magnitude of the star’s magnetic field. This has led 
to the calculation of several theoretical “death lines” in 
the P — P diagram, after crossing which pulsars become 
radio-quiet 8 . These are generally well approximated by 
the equation 

— = 0.17 x 10 12 G s“ 2 (22) 

(jBhattacharva et al.l Il99it) . This theoretical prediction 
is well observed empirically, as there is in fact a well de¬ 
fined cut in the distribution of observed pulsars in the 
P — P diagram (see Figure [4|. In our simulations, equa¬ 
tion [22] will be used to classify synthetic pulsars as ei¬ 
ther active or extinguished radio sources. Actual radio 
pulsar surveys, though, generally have reduced sensitiv¬ 
ity to long-period pulsars, due to radio interference and 
hardware and software high-pass filters. These effects, 
which we do not model, may play a role in the paucity 
of pulsars observed in the long-period part of the P — P 
diagram. We consider in section 14.4.21 whether a sudden 
radio emission shut-off is really required by the data. 

3.4. Radiation Beaming 

The pulsed flux observed from the pulsars is due to 
their highly beamed radio emission, which sweeps our 
line of sight as the star rotates. A consequence of this 
beaming is that only a small fraction of the pulsars 
emit radio w aves in our direction and are observable. 
iTauris fc Manchester! (|l998f ) have analysed polarization 
data for a large number of isolated radio pulsars and de¬ 
veloped a model (TM98) for the fraction / of pulsars that 
are beamed toward us as a function of the spin period: 

f(P) = 0.09[log (P/s) - if + 0.03. (23) 

8 We define the radio-loud pulsars as those that are active radio 
emitters, but that are not necessarily beamed toward us (“poten¬ 
tially observable”). The pulsars that are not radio-loud are said to 
be radio-quiet. 


TABLE 5 

Parkes and Swinburne Multibeam Survey Parameters 


Parameter 


Parkes a Swinburne 8 


Longitude Coverage 
Latitude Coverage 
Aperture Diameter 
Central Observing Frequency ( v) 
Bandwidth (Au) 

Number of Channels (N c h) 
Antenna Gain (G) 
Number of Polarizations (N p ) 
Integration Time (/.,>,/) 
Sampling Interval ( t aa m P ) 
S/N Detection Threshold (rr) 
System Losses Factor (/3) 
Receiver Temperature (Tree) 


260° < l < 50° 

6 <5° 5° < |fc| < 15° 

64 m 
1.4 GHz 
288 MHz 
96 

0.64 K Jy- 1 0 
2 

35 min 265 s 

250 ps 125 ps 

8 

1.5 
21 K 


a Parameters ta ken from [Man chester et alJ 1120011) . b Parame¬ 
ters taken from lEdwards et al.l 11200 11 1. c Average over the 13 
beams of the Parkes Multibeam receiver. 

In this model, the average beaming fraction is ~ 10% 
(both for the pulsar sample analysed by Tauris & Manch¬ 
ester 1998 and for the optimal population model pre¬ 
sented in this work; see section f3.81) . We implement it in 
our simulations by assigning to each Monte Carlo pulsar 
a probability f(P) that its radio beam crosses our line 
of sight. 


3.5. Radio Surveys 

The goal of our population synthesis is to reproduce 
the actual Galactic pulsar population. To assess its suc¬ 
cess, it is necessary to compare the results of our sim¬ 
ulations with observational data. Because the observed 
pulsar sample is heavily biased with respec t to the u n der¬ 
lying population (for a discussion, see e.g. lLorimer et ahl 
[l993h . care must be taken to ensure that equivalent real 
and simulated subsets are compared. We do so by re¬ 
stricting the observed real sample to pulsars detected 
in surveys with well defined sensitivity and performing 
simulations of the same surveys on our synthetic Galaxy. 
We first motivate and precisely define our chosen com¬ 
parison sample, then describe how the radio surveys are 
modelled. 


3.5.1. Comparison Sample 

The ATNF Pulsar Database is currently the most 
comprehensive pulsar catalog available. It contains 
over 1500 objects detected in more than a dozen 
surveys. About 2/3 of the pulsars were detected in 
the PM or the SM surveys at 1.4 GHz. These two 
surveys together covered the well defined sky rectangle 
260° < l < 50°, | 6 | < 15° and used the same observin g 
system ([Manchester^ et al.1 120011 : lEdwards et al.l l2001h . 
Apart from complementary sky coverages (\b\ < 5° 
and 5° < \b\ < 15°) and different integration times 
(35 min vs. 265 s) and sampling intervals (250 /.is vs. 
125 /is), the observing setups are in fact identical, 
providing a remarkably homogeneous pair of surveys 
dominating the available statistics. The parameters 
of the PM and SM surveys are given in Ta ble [5j 
Most o f the earlier su r veys (e. g.. lLarge fc Vaughanl 
19711: iHulse fc Tavloil 119751: lPavies et alJ 1977| : 


Manchester et al 


1978: 


Damashek et all 


Dewev et al.l 119851 : IStokes et al.l I1986T Nice et alJ 
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Parkes Multibecim Survey 
o 



Swinburne Multibeam Survey 
o 



Period (s) Period (s') 

Fig. 5.— I Dewey et alj fll984T ) (dotted) and Crawford’s (dashed) S m in curves for the Parkes (left) and Swinburne (right) Multibeam 
surveys. In each case, the curves are drawn for DM = 0,100, 300,1000 pc cm -3 (in increasing order of Smin). The dots represent pulsars 
detected in the actual surveys that were retained in our comparison sample. The size of the dot indicates whether 0 < DM < 100, 
100 < DM < 300, 300 < DM < 1000, or DM > 1000 pc cm -3 , with the larger dots corresponding to the larger DM. In this figure, Smin 
has been averaged over the 13 beams of the Parkes Mu ltibeam receiver, accoun ting for a Gaussian beam power pattern, sky background 
brightness temperature, and interstellar scattering as in I Manchester et al.l (|20011 ). 


ISaver et~ahlll997t iManchester et al.lll996f ) had a central 
observ ing frequency near 400 MHz. iLorimer et al.l 
(jl995l) have suggested that pulsar spectral indexes may 
be correlated wi th other pr oper ties, in particular the 
spin period. iMaron et al.l (|2000l) , basing their analysi s 
on an extension of the data set of lLorimer et al.l (l995), 
on the other hand, have not found such a correlation. 
In the early stages of this project, we experimented 
with simulations of both 400 MHz and 1.4 GHz surveys, 
assuming spectral indexes independent of other pulsar 
characteristics, and with a normal distr ibution with 
mean —1.8 and standard deviation 0.2 (IMaron et al.l 
l2000f ). Discrepant results were obtained depending 
on whether the pulsar birthrate was calibrated based 
on the 400 MHz or 1.4 GHz detections, suggesting an 
inadequate modelling of the spectral dependence. The 
correct way of simulating spectral variations is unclear. 
Another difficulty in modelling a heterogeneous collec¬ 
tion of surveys is ensuring that the relative sensitivities 
of the surveys are accurately represented. An underes¬ 
timated flux detection threshold for a survey covering 
a given portion of the sky, for example, would result 
in an artificially high number of simulated detections 
in that area. Such spurious features would be difficult 
to distinguish from genuine ones. For these reasons, 
we limit ourselves to the pulsars detected in the PM 
and SM surveys. Of this reduced sample, we further 
ignore the pulsars that lie outside the documented 
survey boundaries, those with P < 30 ms or P < 0, and 
those in binary systems. The latter restrictions limit 
the contamination of our sample by recycled objects. 
The resulting sample contains 1065 pulsars, the PM 


and SM surveys contributing 914 9 and 151, respectively. 
A significant fraction of the retained pulsars (126 and 
97 for the PM and SM surveys, respectively) have 
missing flux or period derivative measurements. The 
pulsars with missing data must nonetheless be taken 
into account to ensure that the PM and SM detections 
are represented in accurate proportion and in order to 
determine the correct number of Monte Carlo pulsars 
necessary to reproduce the number of actual survey 
detections. When comparing the simulated and real 
samples (see section 13.71) . we thus consider all the 
real pulsars, except for the distributions that require 
knowledge of missing data (flux densities, magnetic 
fields, and P — P diagram), for which we assume that 
the subsample with data is unbiased with respect to 
the detections. In the case of the PM survey, there is 
in fact no expected bias (R. N. Manchester and D. R. 
Lorimer 2005, private communication). Unless the bias 
is very important in the case of the SM survey, our 
results are most likely unaffected, as the total subsample 
with missing data represents only 9% of the comparison 
sample. 

3.5.2. Sensitivity Modelling 

The detectability of a pulsar depends on its intrinsic 
properties (brightness, pulse period, and duty cycle), on 
its location (distance, DM, interstellar scattering, and 
brightness temperature of the background sky), and on 
the details of the observing system. The minimum flux 
theoretically detectable from a radio pulsar is usually 

9 Data for 117 PM detections were provided to us prior to public 
release by D. R. Lorimer. 
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estimated by the formula 


, _ ^ P<j(T rec + T s k y ) I We 

1 min — Obeam “ , = = . \ / ~/F. T .. • 

Gy/ N p Avt int V P -We 


(24) 


where Sb eam is a factor accounting for the reduction in 
sensitivity to pulsars located away from the center of the 
telescope beam, T rec is the receiver temperature on cold 
sky, T s k y is the sky background temperature, G is the 
antenna gain, N p is the number of polarizations, Au is 
the receiver bandwidth, t,; nt is the integration time, P 
is the pulse period, W e is the effective pulse width, a is 
the signal-to-noise detection threshold, and 0 is a con¬ 
stant accounting for various system losses dDewev et al.l 
I1984H . Further detail s of the Smin ca lculation are pro¬ 
vided in Appendix 1. ICrawfordl (12000H performed a de¬ 
tailed analysis of the observing system used in the PM 
and SM surveys, providing an alternative, ostensibly 
more accurate, S m in for these surveys. Fi gure [5] shows 
sensit ivity curv es calculat e d usin g both the lDewev et al.l 
(|1984f) and the ICrawfordl (|2000D S m in formulae, for the 
PM and SM surveys. 148 pulsars detected in the ac¬ 
tual PM survey (out of the 806 with flux measurements 
retained in our analysis) have measured 1.4 GHz flux 
< 0.22 mjy, the average nominal survey sensitivity for 
non- millisecond pulsars accor ding to Crawford’s formal¬ 
ism dManchester et al]l2001l f. Similarly, 8 pulsars de¬ 
tected in the SM survey (out of the 54 with flux measure¬ 
ments retained) have measured flux below the average 
survey sensitivity given by Crawford’s S m i n , 0.63 mjy. 
While interstellar scintillation is likely responsible for 
some of these detections, it appears (barring a system¬ 
atic error in the catalogued pulsar fluxes) that the flux 
threshold is overestimated b y Crawfor d’s analy sis. The 
threshold predicted by the iDewev et al.l (Il984tl formula 
is lower and empi rically seems more acc urate. It is thus 
the one we adopt. iGonthier et alJ (j2004l ) used the Craw¬ 
ford formalism in their analysis of the PM data. They 
remarked that t he m o delled sensitiv ity differs from that 
predicted by the lDewev et al.l (jl984h formula, but appar¬ 
ently did not consider which best empirically accounts for 
the actual detections. 

A synthetic pulsar is detected if and only if it lies in 
the sky area covered by either the PM or SM survey and 
the radio flux density of the pulsar at 1.4 GHz at the 
Earth (S 1400 ) exceeds the survey threshold, S m i n . 


3.6. Simulation Procedure 

After describing the components of our population syn¬ 
thesis, we now explain how we proceed with the simu¬ 
lations. We first choose values for the model parame¬ 
ters. Synthetic pulsars are then created one at a time. 
The age (t) of the pulsar is chosen randomly in [0, t m ax\, 
where t max is such that practically all pulsars allowed 
by the model parameters would cross the death line be¬ 
fore reaching that age. More precisely, the spin-down 
equation [20] and the death line equation [22] determine 
the age of radio emission cessation tdeath(Po, B) for a 
pulsar with initial spin period Po and magnetic field B. 
We maximize tdeath over the parameter space rectangle 
[(Po) ~ 3(Tp 0 , (P 0 ) + 3 <7 Po ] X [(logP) - 3CTi ogB , (logP) + 
30iogs] and set t max equal to the maximum. Because 
Po and B are chosen independently, only a negligible 
fraction ~ 0.006 of the synthetic pulsars have a lifetime 


exceeding t max . Typically, t max < 1 Gyr. Astrophysi- 
cally, this is equivalent to assuming a constant birthrate 
in the Galaxy over the life span of the longest lived ob¬ 
jects allowed by the model, i.e. that the population is 
in a steady state. This is a well motivated assumption, 
given the Galaxy’s age of ~ 10 Gyr. Birth characteristics 
are assigned to the pulsar following section 13.21 and are 
evolved as in section 13.31 If the resulting evolved pulsar 
lies beyond the death line or is not beamed toward us, 
we subsequently ignore it (except for keeping count of 
the number of generated MC pulsars, Nmc )• Otherwise, 
we test it for detection in the PM and SM surveys as in 
section 13.5.21 We repeat the procedure until the num¬ 
ber of detections in the simulation equals the number of 
detections in the actual surveys. The estimated pulsar 
birthrate is then N = Nmc/ tmax- The obtained sam¬ 
ple of observed synthetic pulsars is compared with the 
real observed sample to assess the realism of the model 
population. 

3.7. Model Assessment and Search for “Optimal” 
Parameters 

Before adopting a set of model parameters and drawing 
conclusions from it, we must verify that it reproduces the 
actual observations to a satisfactory level. We do so as 
follows. 

For each simulation, the observed marginal distribu¬ 
tions of l, 5, DM , 5i4oo, P, and B in the simulation are 
compared with the equivalent distributions for the real 
sample. We also consider the simulated P — P diagram. 
Except for the surface magnetic field, we have selected 
directly observable quantities that characterize the spa¬ 
tial, rotational, and brightness properties of the pulsars. 

Unfortunately, the comparison does not lend itself well 
to a rigorous, fully quantitative analysis. We have consid¬ 
ered two statistical methods for quantifying our analysis. 

The Kolmogorov-Smirnov goodness-of-fit test, that we 
have already used in section [2j is one of the most pop¬ 
ular tests for comparing statistical distributions. Given 
two empirical cumulative distribution functions F± and 
F 2 for real-valued random variables X\ and X 2 , the 
“D-statistic” is defined as D = sup,„ 6R |Pi(:r) — P 2 (a:)|. 
Remarkably, in the case of the null hypothesis that 
the samples 1 and 2 are identically distributed, the 
.D-statistic has a distribution which is independent of 
that of the compared random variables. We may thus 
test the hypothesis that X\ and A' 2 are identically dis¬ 
tributed by computing the probability of observing a 
deviation D greater than or equal to the one observed 
(D 0 bs), he. the P-value -Pp^-S = -P(D ^ D a b s )- Accord¬ 
ing to the K-S test, we reject the null hypothesis with 
significance level a if Pr-S < a • Two problems prevent 
us from relying exclusively on the K-S test. First, a pul¬ 
sar is not adequately characterized by a single real num¬ 
ber as its location, spin characteristics, and brightness 
are all of interest. It is a mathematical fact that the K-S 
D-statistic loses its distribution-free virtue for multidi¬ 
mensional random variables and no p ractical generaliza - 
tion to higher dimensions is available (lJustel et al.lll997l ). 
Alternatively, we may perform a K-S test on each of the 
marginal distributions of interest. Because they are not 
independent, we may not however meaningfully combine 
the results to obtain a global goodness-of-fit figure for 
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the model. The second problem with the K-S test, in 
our context, is that it tests whether the two samples are 
drawn from exactly the same distribution. In the case 
of large samples, the test is sensitive to small deviations 
and, if present, they often lead to a formal rejection of 
the null hypothesis to a very low significance level. Such 
deviations may turn out to be of little astrophysical im¬ 
portance and are expected given the complexity of the 
system that we are attempting to model in a relatively 
simple way. Some human judgement is thus required 
when interpreting the quantitative results. 

Bayesian maximum likelihood ana lyses, such as the one 
performed by lArzoumanian et ~akl ( 20021 or our treat¬ 
ment of the birth velocity distribution in section[2J do of¬ 
fer a procedure for optimizing a given parameteric model 
and are not limited to unidimensional distributions. A 
major disadvantage, though, is that they do not provide 
any absolute measure of the correctness of the assumed 
functional form. Additional checks, e.g. visual and using 
K-S tests as we did in section [2j are required for that 
purpose. Because of this, we choose not to tackle the 
computational challenges associated with the implemen¬ 
tation of maximum likelihood formalism for multidimen¬ 
sional distributions. 

In light of these remarks, we adopt a semi-quantitative 
approach. We start with artibrary, but reasonable, 
model parameters and iterate. At each step, we perform 
independent K-S tests on the marginal distributions (of 
l, b , DM. 5*1400, P, and B ) and use the K-S P-values to 
guide our search for “optimal” model parameters. The 
choice of subsequent model parameters, though, is based 
on physical intuition and involves qualitative considera¬ 
tions. We stop when no significant further improvement 
appears possible. We do not have a well defined metric 
for the “best” model parameters and there is no guar¬ 
antee that we have in fact found them. We do however 
believe that substantial improvement would require revi¬ 
sion of the assumed functional forms in the population 
model. 

Hereafter, we commit the abuse of language of referring 
to the final model as “optimal” or “best”, cautioning that 
we do so very informally. 

3.8. Results 

Tabic [HI lists the parameters of the optimal model. The 
marginal distributions of observed l, b , DM, Shoo, P, 
and B for the model are shown in Figured] For the sim¬ 
ulation histograms, the number of pulsars in each bin is 
the average over 50 MC realizations of the model. The er¬ 
ror bars indicate the corresponding standard deviations. 
On each histogram, the associated K-S P-value is dis¬ 
played in the upper right corner. For its computation, 
all 50 x 1065 MC pulsars are used to construct the empiri¬ 
cal distribution function for the model. A MC realization 
of the P — P diagram for the model is shown in Figure 0 
The agreement with the actual distributions and P — P 
diagram (c.f. Figure [4} is generally qualitatively good. 
We discuss the results of our population synthesis and 
their implications in the next section. 

4. DISCUSSION 

4.1. Spatial Distribution of the Pulsars 

We have taken the first step in modelling the rich 
Galactic structure by simulating the four major spiral 


TABLE 6 

Optimal Population Model Parameters 


Model Parameter 

Value 

Radial Distribution Model 

Yusifov &; Kugiik 

Tii 

0.55 kpc 

a 

1.64 

b 

4.01 

Birth Height Distribution 

Exponential 

{^0> 

50 pc 

Birth Velocity Distribution 

Exponential 

(V3D) 

380 km s” 1 

Birth Spin Period Distribution 

Normal 

(Po) 

300 ms 

°Pn 

150 ms 

Magnetic Field Distribution 

Log-Normal 

(log (B/G)> 

12.65 

^log B 

0.55 

Luminosity Model 

P — P Power Paw 

Lo 

0.18 mjy kpc 2 

ep 

-1.5 

ep 

0.5 

Cr L CO rr 

0.8 


arms and a non-trivial Galactocentric radial distribution 
of pulsar birth sites. 

4.1.1. Galactocentric Radial Distribution 

Following the completion of the first high-frequency 
(1.4 GHz) ra dio pulsar surveys capable of probing th e 
inner Galaxy dcTifton et al.lll992t I Johnston et al.lfl992f) . 
iJohnstoiil (119941) concluded that the Gaussian radial dis¬ 
tribution models peaking at the GC often assumed in 
pulsar population studies were incompatible with the ob¬ 
servations. Instead, he proposed a model in which the 
surface density of pulsars peaks at a distance of 4 kpc 
from the GC. In spite of this, the practice of assuming 
simplistic radial distributions in pulsar population stud¬ 
ies has generally persisted. 

Analy s ing re cent data from the PM and S M surveys, 
iLorimed (|2003l) and lYusifov fe Kiiciikl (|2004D reiterated 
the pulsar deficit near the GC. In our population syn¬ 
thesis, we have ad opted the Galacto c entric radial distri¬ 
bution derived bv lYusifov fe Kiichkl (|2004T ) for the true, 
evolved Galactic population as the distribution for the 
birth locations of our synthetic pulsars. We first exam¬ 
ine whether it is indeed consistent with the observations. 

The observed distributions of pulsars in Galactic lon¬ 
gitudes (as angular position indicators) and in dispersion 
measures (as distance indicators) provide a test for the 
assumed radial distribution of the synthetic pulsars. The 
agreement between the simulated and real distributions 
is reasonably good in both cases (see Figure [6]), the ma¬ 
jor features (fall-off in detection density for \l\ > 40° 
and approximate plateau in-between; rapid rise of the 
DM distribution for DM < 200 pc cm” 3 and slow de¬ 
cline afterward) being well reproduced. In the case of the 
Galactic longitudes, a noticeable discrepancy is the lack 
of a significant decrease in simulated detection density for 
|/| < 20° that is observed in reality and which suggests 
spiral arm tangents near l = ±30°. The source of the dis¬ 
crepancy is unclear. Three possibilities are an inadequate 
modelling of the spiral arm structure (see section 14. 1.21) , 
an incorrectly modelled Galactocentric radial distribu¬ 
tion, or underestimated selection effects against pulsars 
near the GC (e.g., underestimated scattering). Unfortu- 
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Galactic Longitude (°) 


Galactic Latitude (°) 





Period (s) Surface Magnetic Field (G) 

Fig. 6. — Distributions of observed pulsar Galactic longitudes and latitudes, dispersion measures, flux densities at 1.4 GHz, pulse periods, 
and surface magnetic fields for our optimal model (solid lines) compared to the real distributions (hatched histograms). For the simulation 
histograms, the number of pulsars in each bin is the average over 50 Monte Carlo realizations of the model. The error bars indicate 
the corresponding standard deviations. On each histogram, the associated K-S P-value is displayed in the upper right corner. For its 
computation, all 50 X 1065 MC pulsars are used to construct the empirical distribution function for the model. 


nately, we cannot isolate or exclude any of them on the 
basis of our analysis, for the parameters of our model are 
too intertwined with each other to single out the effect 
of each one. In the case of the dispersion measures, the 
simulated distribution appears somewhat skewed toward 
higher values with respect to the real o ne. This is in 
qualita tive agreement with a critique by iKramer et al.l 
(|2003aj) of the NE2001 model. They considered the dis¬ 
tribution of derived distances from the Galactic plane 
(z) of known pulsars with \b\ < 20° as a function of dis¬ 
tance from the Sun using NE2001. They found a slight 
decreasing trend, which is opposite to what is expected, 
as the maximum 2 -distance probed increases with dis¬ 
tance from the Sun. This suggests that NE2001 has 
a tendency to underestimate the distances to pulsars, 
i.e. to overestimate the free electron density. The minor 


discrepancies between the real and simulated distribu¬ 
tions appear equally possible as being due to remain¬ 
ing imperfections in the NE2001 model rather than the 
radial distribution model. Thus, within the limitations 
imposed by our method, the Galactocentric radial distri- 
bution of pulsar b i rth sit es is consistent with the model of 
lYusifov fe Kiiciikl (l2004f) and has a deficit near the GC. 

We now addr ess the question o f whet her it was jus¬ 
tified to use the lYusifov fc Kiiciikl (120041 ) model for the 
evolved pulsars as the birth distribution in the first place. 
In Appendix 2, we compute the birth radial distribu¬ 
tion, taking the lYusifov fc Kiiciikl (1 2004 ) model for the 
evolved distribution as a prescription, under the assump¬ 
tion that the kinematic and rotational evolution of the 
pulsars is accurately described by our optimal model pa¬ 
rameters. The top panel of Figure [8] shows the observed 
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Period (s) 

Fig. 7.— P — P diagram for a typical Monte Carlo realiza¬ 
tion of our optimal model. Each point represents a simulated pul¬ 
sar. Contours of constant surface magnetic field and of constant 
characteristic age are indicated by dashed and dotted-dashed lines, 
respectively. The thick solid line marks the modelled death line. 


distribution in Galactic longitudes obtained in simula¬ 
tions when this birth distribution is assumed and all the 
other parameters are as in the optimal model. There 
are clear peaks in the detected simulated pulsar den¬ 
sity at Galactic longitudes ~ 295°, ~ 215°, and ~ 30°. 
Moreover, the general “U” shape of the distribution is 
strongly suggestive that pulsars are more preferentially 
detected in a thin ring around the GC in the simula¬ 
tions than in reality. The most plausible explanation 
is the that birth radial distribution assumed (see Fig¬ 
ure [9| is too strongly concentrated in such a ring. While 
this suggests an in consisten c y betw een our optimal model 
and th e studies of lLorimeil (1200311 and lYusifov fc Kiiciikl 
(i200l . we argue that this discrepancy may easily be un¬ 
derstood if the radio luminosities of pulsars are a sim¬ 
ple function of their period and period derivative such 
that younger pulsars have a tendency to be more lumi¬ 
nous and therefor e preferentially detec ted. In fact, if this 
is the case, then iLorimerl ([2003) and lYusifov fe Kiiciikl 
(|2004T) . who have not treated this bias in their studies, 


may have derived radial distributions that are biased to¬ 
ward younger pulsars. If this bias is important, then their 
distributions may actually be better approximations to 
the birth radial distribution than to the unbiased evolved 
one. In any case, the survey results are more closely re- 
pro duced in the simu l ations presented here when taking 
the lYusifov fe Kiiciikl (I2004H distribution as the birth ra¬ 
dial distribution than with the more concentrated dis¬ 
tribut ion inferred from prescribing the Yusifo v~ fe Kiiciikl 
(12004) model for the evolved pulsars, and this motivates 
us to adopt the former case. We shall see in section l4~4l 
that the data strongly suggest that the radio luminosity 
of pulsars must indeed be somehow correlated with their 
age. If this is correct, it may also explain discrepancies 
regarding the torque decay in isolated pulsars between 
various studies in the literature (see section H31) . 

To further verify that the consistency of the 


Birth Radial Distribution Estimated from the 
Evolved Model of Yusifov & Kucuk (2004) 



Gaussian Birth Radial Distribution 



Galactic Longitude (°) 


Uniform Birth Radial Distribution 



Fig. 8 .— Distributions of observed pulsar Galactic longi¬ 
tudes (solid lines) compared to the real distributions (hatched his¬ 
tograms) for alternative models of the birth radial distribution of 
pulsars. The top panel shows the distribution obtained when the 
lYusifov Sz Kiiciikl (120041) model is taken as a prescription for the 
evolved distribution and the corresponding birth distribution is es¬ 
timated as in Appendix 2. Clear peaks in the detected simulated 
pulsar density are seen near ~ 295°, ~ 215°, and ~ 30° and the 
general “U” shape of the distribution is suggestive that pulsars 
are more preferentially detected in a thin ring around the Galactic 
Center in the simulations than in reality. In the center panel, a 
Gaussian distribution peaking at the GC and with scale radius 5 
kpc was assumed. An excess of simulated detections is seen near 
the GC (|Z| < 20°). In the bottom panel, an uniform distribution 
on a disc of radius 15 kpc centered on the GC was assumed. An 
excess of simulated detections is seen for |/| > 40° and there is 
a paucity elsewhere. In each case, all other parameters are as in 
the optimal model. For the simulation histograms, the number of 
pulsars in each bin is the average over 50 Monte Carlo realizations 
of the model. The error bars indicate the corresponding standard 
deviations. On each histogram, the associated K-S P-value is dis¬ 
played in the upper right corner. For its computation, all 50 x 1065 
MC pulsars are used to construct the empirical distribution func¬ 
tion for the model. 


lYusifov fe Kiiciikl (I2004T) model used as the birth radial 
distribution of pulsars is not an artefact of the poten- 
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Distance from the Galactic Center (kpc) 

Fig. 9. — Probability density function fo r the evolved Galactocen- 
tric radial coordinate of a pulsar for the lYusifov Sz Kiiciild (|2004D 
model (solid line) and estimate of the corresponding birth PDF 
as computed in Appendix 2 (dashed line). The shaded histogram 
shows the empirical evolved distribution for simulated pulsars with 
birth radial distribution given by the dashed curve. All the other 
simulation parameters are as in the optimal model. 

tial insensitivity of our simulation results to the details 
of the birth spatial distribution of the pulsars and that 
the complexity of the model is justified, we have con¬ 
sidered two simpler models. Namely, we have tried a 
birth radial distribution in which the pulsar surface den¬ 
sity decreases like a Gaussian with scale radius 5 kpc 
with distance from the GC and one in which the pulsar 
density is uniform on a disc of radius 15 kpc centered 
at the GC, keeping the rest of our optimal simulation 
parameters. The resulting distributions in Galactic lon¬ 
gitudes are shown in the center and bottom panels of 
Figure [8] A clear excess of simulated detections near the 
GC (|6| < 20°) is seen in the Gaussian case. In the uni¬ 
form case, there is an excess away from the GC (b > 40°) 
and a paucity elsewhere. In both cases, these discrepan¬ 
cies are qualitatively expe cted from a compa r ison o f the 
assumed models with the lYusifov fe Kiiciild (120041 1 dis¬ 
tribution and so the complexity of the latter model well 
motivated. 

As discussed bv lYusifov fc Kiictikl (j2004h . the modelled 
radius of maximum pulsar density, 3.15 kpc, is smaller 
by nearly 1.5 kpc than the peak r adii o f several Popu¬ 
lation I objects. iBronfman et alJ (|2000i) found 4.7 kpc 
and 4.3 kpc for the radii of maximum far infrared lumi¬ 
nosity produced by massive stars and of the H 2 surface 
density, respectively. The s h ell supernova remna nt dis¬ 
tribution inferred bv lGase fc Bhattacharval (Il998t ) peaks 
at a distance of about 4.8 kpc from the GC. Nevertheless, 
the qualitative concentration in a ring midway between 
the Sun and the GC is common to both the pulsars and 
Population I objects. The quantitative discrepancy may 
again be attributable to remaining imperfections in the 
modelling of the free electron density in NE2001. How¬ 
ever, if NE2001 indeed has a tendency to underestimate 
distances, then the distance to the pulsar ring most likely 
has been underestimated, exacerbating the discrepancy 
with the location of massive stars. While the origin 
the quantitative discrepancy is unclear, the uncertain¬ 
ties remaining in the determination of the exact pulsar 
radial distribution prevent us from concluding that the 
disagreement is genuine. As the birth sites of pulsars 
are also correlated with the spiral arms (section 14.1.21) , 
active star-forming regions, and their scale height is also 


r c <2 Myr 



Galactic Longitude (°) 
r c >2 Myr 



Fig. 10. — Distributions of observed pulsar Galactic longitudes 
for our optimal model (solid lines) compared to the real distribu¬ 
tions (hatched histograms), for the pulsars with characteristic age 
r c < 2 Myr (top) and for those with r c > 2 Myr (bottom). For the 
simulation histograms, the number of pulsars in each bin is the 
average over 50 Monte Carlo realizations of the model. The error 
bars indicate the corresponding standard deviations. On each his¬ 
togram, the associated K-S P-value is displayed in the upper right 
corner. For its computation, all 50 x 1065 MC pulsars are used to 
construct the empirical distribution function for the model. Near 
±30° (corresponding to the Norma and Perseus, and to the Crux- 
Scutum and Carina-Sagittarius arms, respectively), there appears 
to be an enhanced pulsar detection density in both the real and 
simulated samples in the sample of young pulsars, although the 
simulated distribution lacks a pronounced dip near the Galactic 
Center. For the pulsars older than 2 Myr, the correlation between 
the spiral arms and the Galactic longitudes is greatly diminished. 

consistent with that of massive stars (see 14.1.31) . it is safe 
to assume that, as theoretically expected, pulsar birth 
locations do coincide with those of massive star deaths. 

4.1.2. Spiral Arm Structure 

From the Galactic longitude distribution of our optimal 
model in Figure dl on which we have indicated the tan¬ 
gent points to the modelled spiral arms (c.f. Figured]), 
we see that near ±30° (corresponding to the Norma and 
Perseus, and to the Crux-Scutum and Carina-Sagittarius 
arms, respectively), there is a hint of enhanced pulsar 
detection density in both the real and simulated sam¬ 
ples. This correlation is more apparent when only young 
pulsars (r c < 2 Myr) are considered (see Figure [TOl). al¬ 
though the simulated distribution lacks a pronounced dip 
near the GC, as in section 14.1.11 For the pulsars older 
than 2 Myr, the correlation between the spiral arms and 
the Galactic longitudes is greatly diminished. 

Interestingly, the correlation between the young simu- 
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lated pulsars and the tangent points near —50° and —72° 
seems stronger than between the real young pulsars and 
the tangents. This may be due to an important assump¬ 
tion that we have made in modelling the spiral arms. 
We have assumed in our simulations that the spiral pat¬ 
tern is fixed in time in the coordinate system corotat¬ 
ing with the Sun. In reality, the matter in the Galactic 
disc does not in general co r otate with the spiral pat¬ 
tern (see, e.g.. lBinnev fe TremaineIl987f) . In our Galaxy, 
the corotation radius is ~ 14 kpc f Burtoniri97lh . so that 
from the vantage point of the Sun, the spiral pattern 
and hence the presumed birth sites of pulsars move with 
time. Thi s effect is significa nt , as w as demonstrated by 
iRamachandran fe Deshpandel (|1994f) who found a corre¬ 
lation between the current pulsar distribution and the 
expected locations of the spiral arms ~60 Myr ago, con¬ 
firming that the relative angular velocity between the 
spiral pattern and the Sun motion spreads the pulsars 
over the disc. That this effect was not included in the 
simulations may explain the stronger correlation between 
the current spiral arm tangents and the young simulated 
pulsars. That the correlation between the real young pul¬ 
sars and the spiral arms near l = ±30° is preserved may 
due to the fact that the Norma and Perseus, and the 
Crux-Scutum and Carina-Sagittarium arms are tightly 
wound in those regions, acting as larger concentrations 
that take longer to move a distance exceeding their size. 

Thus, while spiral arm structure is certainly required 
to reproduce the spatial distribution of pulsars, our sim¬ 
ple modelling using fixed spirals appears deficient. A 
treatment of the kinetics of the spiral pattern may be 
needed for more accurate results. 

4.1.3. Scale Height 

The distribution of distances of the pulsars from the 
Galactic plane is most directly reflected in the observed 
distribution in Galactic latitudes (Figure [6]). We find 
that (zq) = 50 pc, corresponding to a e _1 distribution 
thickness of 100 pc, accomodates the observations well. 
However, due to the high birth space velocities of the 
pulsars, our simulations are not sensitive to the mean 
birth distance from the Galactic plane to a precision bet¬ 
ter than a few tens of parsecs. Massive star formation 
is distributed on a disc l ayer of thickness ~70 pc (full 
width at half maximum; iBronfman et al.l 12000 1. Using 
a vertical velocity 10 km s^ 1 and a main sequence life¬ 
time 10 Myr characteristic of massive stars, these move 
~ 100 pc away from their formation site perpendicular 
to the Galactic plane before exploding in a supernova, 
somewhat thickening the expected spatial distribution of 
these events, but not enough to create any serious dis¬ 
crepancy with the birth scale height of pulsars in our 
simulations. Thus, up to the precision provided by our 
simulations, the birth scale height of pulsars coincides 
with that of massive stars ending their lives in super¬ 
novae. 

4.2. Birth Spin Period Distribution 

The birth spin period of the Crab pulsar, whose age 
is known from its association with the supernova of 
A.D. 1054 and whose braking index has been measured, 
was th e first to be est i mate d, with a value Pq ^19 ms 
(iManchester fe Tavloilll977lh This has led to the gen¬ 
erally made assumption that the initial spin periods of 


pulsars are much smaller than observed ones. The exis¬ 
tence of a young pulsar with period 16 ms has p rovided 
additional support to that assumption (Marshall ct al. 
Il998ll2004h . However, a number of recent measurements 
suggest that initial periods in the range ~ 50 — 150 ms 
are not uncommon (see Table 0 which lists the pulsars 
with estimated birth spin periods). In addition, it has 
been argued that the 105 ms pulsar PSR J1852—0040 
may have had a birt h period close to the present value 
dGotthelf et al.ll2005ll . If, as our results suggest (see sec¬ 
tion EH]), short-period pulsars tend to be brighter, it is 
possible that only the lower part of the birth spin period 
distribution has been sampled and it may extend to con¬ 
siderably larger values. Several population studies have 
supported this idea, either assuming or requiring a large 


fraction of pulsars to be born with periods > 100 ms 

fe.g.. IVivekanand & Naravanl 119811: Naravan! 

1987c 

Emmering & Chevalier 1984 iNaravan & Ostriker 

”991 

Regimbau & de Freitas Pacheco 20011; Gonthier et al. 

2004|). From a pulsar current analysis of the PM 

data, 

many 

Vranesevic et al.l (2004T) estimated that perhaps as 


as 40% of the pulsars may be born with periods in the 
range 0.1 — 0.5 s. In our optimal model, the birth spin 
period distribution is normal, centered at 300 ms and 
with standard deviation 150 ms. Unfortunately, this dis¬ 
tribution is not precisely constrained by our method, as 
the spin period of a pulsar loses the imprint of its ini¬ 
tial value as the pulsar ages, being asymptotically de¬ 
termined only by the star’s age and magnetic field (Eq. 
EH). Moreover, the derived distribution is clearly depen¬ 
dent on the details of the assumed luminosity law, as 
we favor one in which the luminosity is dependent the 
spin period of the pulsar (see section lT4l) . For example, 
our simulations are incompatible with a large fraction of 
pulsars born with periods < 100 ms, since these are gen¬ 
erally more luminous and easily detected in our P — P 
luminosity model, resulting in an excess of detections at 
short periods in the simulations with respect to the ac¬ 
tual observations. Nonetheless, as we will argue, the data 
do put constraints on the possible luminosity laws which 
suggest that our results are probably at least qualita¬ 
tively robust. We do not find evidence for a multimodal 
spin period distribution (injection), although we cannot 
r ule it out. 

iTorii et akl (|1999t ) and iKasoi et al.l (|2001h have noted 
that PSR J1811—1925 inside the supernova remnant 
Gil.2—0.3 has an actual age (inferred from a probable 
association with a supernova witnessed by Chinese as¬ 
tronomers in A.D. 386) smaller by a factor ~12 than its 
characteristic age obtained from timing (24,000 yr). An 
important assumption in the determination of the age of 
a pulsar from timing is that its initial spin period is neg¬ 
ligible with respect to the obser ved one ( Pq/P 1; see, 
e.g JTavlor fc Manchesteii[l977t ). A discrepancy between 
the actual and characteristic ages suggests that this as¬ 
sumption is incorrect and that the characteristic age is a 
poor indicator of the pulsar’s true age. As the distribu¬ 
tion of pulsar birth spin periods appears to extend much 
above 100 ms, that characteristic ages are poor age in¬ 
dicators for young pulsars may turn out to be the rule 
rather than the exception. In Figure fill we have plotted 
the fractional age difference At/t rea i = (r c — t rea i)/t rea i, 
where t rea i is the true age of the pulsar, for pulsars with 
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TABLE 7 

Estimated Birth Spin Periods of Young Pulsars 


Pulsar Name Estimated Birth Spin Period References 

(ms) 


J0205+6449 

60 

1 

B0531+21 (Crab) 

19 

2 

J0537—6910 a 

~11 

3, 4 

J0538+2817 

139 

5, 6 

B0540-69 

30 ± 8 

7 

J1124—5916 

>90 

8 

J1811—1925 a 

62 

9, 10 

J1833—1034 

>55 

11 

B1951+32 

27 ± 6 

12 


References. — _ (1) IMurrav et al.l 12002ft : _ (2) 

IManchester^&^a^IoiJ |J1977|) ; ( 3]__|MarshaU^t aL |_j 199811 : 

V4)_JMarsha!ret aL| ll2004l~) : (o) IKra mer et all (j2003h|j i __(6) 
IRomani h Nel < 2003fl (7 ) IRevnoldsI (11 9851): f8~> ICamil o et al.l 
<20021: 191 ITorii et al.l <19991: (101 IKaspi et al.l <20011) : flTT 
ICamilo et al.1 (120051 1 : (121 IMigliaazojetEjJ Jj20021 
a Not observed as a radio pulsar IlCrawford et al.ll 199811 ■ 



treai < 100, 000 yr in a MC realization of our optimal 
model. For given birth spin period Pq and magnetic field 
B , the spin-down equation 1201 implies 


where 


At 


treai 


v{Po,B) 


v(Po,B) 


3Ic 3 Pq 

167r 2 i? 6 R 2 sin 2 x' 


(25) 

(26) 


In particular, the characteristic age with this spin- 
down law is a systematic overestimate of the true age. 
The median value of r](Po,B) for our optimal model is 
~ 60, 000 yr. As the Figure shows, pulsars with true age 
< 30, 000 yr have median characteristic age differing by 
a factor significantly larger than unity from the true age, 
if the optimal model presented here accurately represents 
the true pulsar population. 

We have assumed that the spin down of pulsars is dom¬ 
inated by magnetic braking. Some work has suggested 
that a significant amount of rotational kinetic energy 
in young neutron stars may be carried away by grav- 
itatio nal waves generated by unstable r-mode oscilla- 


tions 

Andersson 1998; Lindb 

om et al. 

1998; Owen et al. 

1998; Andersson et all 19991: 

Ho & Lai 

200®. The spin 


down due to gravitational radiation by a neutron star is, 
however, expected to be important only durin g the first 
year after its formation (|Lindblom et al.lll998l ) and more 
recent investigations suggest that r-modes may saturate 
at low amplitudes a nd therefore not be as important 
as initially thought (|Gressrrian et al.l 120021 : lArras et al.l 
120031: [Brink et al.ll2004j) . Thus, for the purpose of age de¬ 
termination, it appears sufficient to a very good approxi¬ 
mation to consider only the magnetic spin-down history. 
What we have been calling the birth spin period, though, 
is more accurately defined as the period at which mag¬ 
netic braking becomes dominant. 


4.3. Birth Velocity Distribution 

In our population synthesis, we have assumed the ex¬ 
ponential birth velocity distribution favored in section 0 
although it was shown that a two-com ponent Ga ussian 
and single-component Lorentzian and iPaczvnskl (119901 ) 
models were equally consistent with the proper-motion 


TrUe Age (yr) 

Fig. 11.— Fractional difference between the characteristic age 
and the true age of synthetic pulsars younger than 100,000 yr in a 
Monte Carlo realization of the optimal model as a function of true 
age. The solid curve (60,000 yr/true age) represents the model 
median value. 


data considered. To estimate the sensitivity of our re¬ 
sults to this choice, we repeated our “optimal” simula¬ 
tions, but replacing the exponential model with each of 
the other models in turn. The results were by and large 
unaffected. Figure [12] shows the observed distribution 
of Galactic latitudes, which is the most sensitive to the 
assumed birth velocity model, in each case. 

4.4. Radio Luminosity 

4.4.1. Inadequacy of Random Luminosities 

The random luminosity model is clearly inadequate. 
While it reproduces the observed flux distribution rea¬ 
sonably well in comparison with our optimal model (c.f. 
Figure it leads to an observed scale height that is 
too large (Figure fl3l). Moreover, it produces a clear pile- 
up of observed objects near the death line in the P — P 
diagram (Figure fl4l) that is not seen in reality (c.f. Fig¬ 
ure SJ. These discrepancies are largely unaffected by the 
choice of parameter values. These two major discrep¬ 
ancies suggest that, more generally, random luminosities 
lead to the predicted detection of too many old pulsars 
and that the correct model must favor the detection of 
young objects. 

4.4.2. Period-Period Derivative Power Law 

An obvious way of correlating the pulsar luminosity 
with age is to make it a function of the pulsar’s period 
and period derivative. When adjusting the model pa¬ 
rameters for the P — P power-law model, we found that 
in order to simultaneously reproduce the distributions 
in observed periods and magnetic fields and the P — P 
diagram, the power-law exponents ep and ep must be 
close to —1.5 and 0.5, respectively. The simulations did 
not allow us to constrain these parameters to a preci¬ 
sion better than a few tenths and we adopted the round 
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2 —Component Birth Velocity Distribution 



Lorentzian Birth Velocity Distribution 



Galactic Latitude (°) 


Paczynski Birth Velocity Distribution 



Fig. 12.— Distributions of observed pulsar Galactic lati¬ 
tudes (solid lines) compared to the real distributions (hatched his¬ 
togr ams), for Gaussia n two-component (top), Lorentzian (center), 
and IPaczvnskH d 199QD (bottom) birth velocity distributions with 
parameters derived in section [2] All other parameters are as in 
the optimal model. For the simulation histograms, the number of 
pulsars in each bin is the average over 50 Monte Carlo realizations 
of the model. The error bars indicate the corresponding standard 
deviations. On each histogram, the associated K-S P-value is dis¬ 
played in the upper right corner. For its computation, all 50 x 1065 
MC pulsars are used to construct the empirical distribution func¬ 
tion for the model. Each birth velocity distribution is seen to agree 
with the observational data similarly well as the exponential one 
assumed in the optimal model (c.f. Figure |6}. 


numbers. To illustrate the wide range of pulsar lumi¬ 
nosities consistent with the observed sample, we have 
plotted the underlying distribution of luminosities at 1.4 
GHz in a MC realization of our optimal model in Fig¬ 
ure 1151 It is well fit by a log-normal distribution with 



Galactic Latitude (°) 



Flux Density at 1.4 GHz (mdy) 

Fig. 13.— Distributions of observed pulsar Galactic latitudes 
and flux densities at 1.4 GHz for our optimal model with the P — P 
power-law luminosity model replaced by the random model (solid 
lines) compared to the real distributions (hatched histograms). For 
the simulation histograms, the number of pulsars in each bin is 
the average over 50 Monte Carlo realizations of the model. The 
error bars indicate the corresponding standard deviations. On each 
histogram, the associated K-S P-value is displayed in the upper 
right corner. For its computation, all 50 x 1065 MC pulsars are 
used to construct the empirical distribution function for the model. 
The observed Galactic latitude distribution in the simulation is 
clearly thicker than the real one. The sharp drop near |6| =5° is 
due to the reduced sensitivity of the Swinburne Multibeam survey 
relative to the Parkes Multibeam survey. The agreement between 
the flux density distributions is comparable as in the optimal case 
with the P — P power-law model (c.f. Figure [6]). 


mean (logL) = —1.1 (L = 0.07 mJy kpc 2 ) and standard 
deviation ai og L = 0.9 in the logarithm to the base 10. 

We note two interesting facts. First, since B oc y/PP, 
the modelled luminosity contour lines (ignoring the cor¬ 
rections) are parallel to the death line (Eq. 1^1) . Thus, in 
this model, the pulsar luminosity is a function of the dis¬ 
tance from the death line in the P—P plane and they uni¬ 
formly fade down as they approach it. This explains why 
we do not observe a pile-up near the model death line. In 
a simulation with the optimal parameters listed in Table 
[G] but with the death line not applied and t max = 15 Gyr, 
30 ± 5 pulsars are detected beyond the line, compared to 
only 5 in the real sample considered (c.f. Figured]). This 
suggests that a complete radio emission shut-off near the 
modelled death line is nevertheless required, although it 
is inconclusive due to the unmodelled reduction of sensi¬ 
tivity to long-period pulsars (see section [3.3.3D . Second, 
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Period (s) 

Fig. 14. — P — P diagram for a typical Monte Carlo realization 
of our optimal model with the P — P power-law luminosity model 
replaced by the random model. Each point represents a simulated 
pulsar. Contours of constant surface magnetic field and of constant 
characteristic age are indicated by dashed and dotted-dashed lines, 
respectively. The thick solid line marks the modelled death line. 
There is a clear pile-up of pulsars near the death line that is not 
seen in the real diagram (c.f. Figure nil. 



Luminosity at 1 ..4 GHz (mjy kpc 2 ) 

Fig. 15. — Underlying distribution of pulsar luminosities at 
1.4 GHz in a Monte Carlo realization of our optimal model. The 
sample consists of all synthetic pulsars that have not crossed the 
model death line. The solid curve represents the best-fit log-normal 
distribution, with mean (logL) = —1.1 (L = 0.07 mjy kpc 2 ) and 
standard deviation cr\ 0 gL = 0-9 in the logarithm to the base 10. 

the luminosity contour lines are furthermore coincident 

with those of the quantity v7, where E = An 2 IP/P 3 
(the “spin-down lumin osity”) is the rate of loss of ro¬ 
tational kinetic energy ijLorimer fc Kramein2005[ ). This 
in turn is proportional to the potential drop available 
for particle accelera tion in magnetosphere models (e.g., 
iGoldreich fc .Julian! ll969( 1. which may be related to the 
physics of radio emission (see, e.g.. IStollmanlll987bl who 
reached conclusions similar to ours). 

From t heir maximum likel ihoo d analysis of 400 MHz 
surveys. [Arzoumanian et ah! (200 2) inferred a total phys¬ 
ical luminosity also scaling approximately as \/T (more 
precisely, as However, in their 

model, this corresponds to the integrated radiated power 


over Gaussian core and conal beams, whose size and rel¬ 
ative contributions are functions of the pulsar’s period. 
The observed flux, then, is not in general proportional 
to p-i-Spo.^ but depends on t he viewing geometry . 
We co nsidered the modifica t ion bylGonthier et ahl (l200ll ) 
of the I Arzoumanian et~akl (120021) luminosity model, in 
which spectral dependence was included to model the lu¬ 
minosity at 1.4 GHz. For the simplest case in which we 
observe an orthogonal rotator whose beam is confined 
to a plane containing the line of sight, the flux averaged 
over the spin period of the pulsar scales approximately 
as p-0- 8 p0- 4 . The flux contour lines then are no longer 
parallel to the death line and the modelled pulsar fluxes 
are thus quite different than ours. Moreover, the flux 
variations due to beaming geometrie s span several o rders 
of m agnitude (see, e.g., Figure 2 of i Arzoumanian et all 
[2002i) . whereas in our case, the typical correction factor 
due to the imposed dither is ~ 10. 

The simple pseudo-luminosity law that we have used 
may be interpreted as what would effectively result if 
we had geometrically modelled the pulsar beams as be¬ 
ing flat and sharp-edged, so that the observed flux is in¬ 
dependent of the impact angle, provided that the beam 
crosses our line of sight. The dithering would then mostly 
account for intrinsic variations about the standard can¬ 
dle law. iTauris fe Manchester! (119981 ). whose beaming 
model we have adopted, have indeed found in their anal¬ 
ysis of polarization data that the fluxes of pulsars are 
consistent with being independent of t he impact a ngle, 
contra ry to the luminosity models of lArzoumanian et all 
(|2002t ) and iGonthier et al.l (120041 ). This supports the 
hypothesis that emission within the beam boundary is 
patch y, with a random distribution of component loca¬ 
tions dLvne fe Manchesteilll988HManchesteilll995l) . 

4.5. Spin-Down Torque and Magnetic Field 

We have throughout assumed that the magnetic field 
of our synthetic pulsars was dipolar, constant, and ne¬ 
glected possible evolution of the magnetic inclination an¬ 
gle. Since the observed data set seems to be well repro¬ 
duced, our simulations suggest that no significant torque 
decay occurs in isolated pulsars over their lifetime as 
radio-loud sources (~ 100 Myr). In particular, this sug¬ 
gests that magnetic field decay is not significant over this 
ti me scale in isolated ne utro n stars. _ 

IGonthier et al.l (120021 ) and IGonthier et ahl (120041 ) per¬ 
formed simulations similar to ours and, in contrast, found 
evidence in favor of field decay on time scales as short 
as ~ 2.5 — 5 Myr. The main argument, in both cases, 
is that without field decay, a pile-up of simulated de¬ 
tections is seen near the death line in the P — P dia¬ 
gram. Our results suggest that the need for field de¬ 
cay may be an artefact of the luminosity models they 
have assumed. While the first study effectively assumed 
a dithered pseudo-luminosity law like ours and the sec¬ 
ond a seemingly diffe rent , more deta iled geometric model 
based on that of lArzoumanian et all (I2002D . in both cases 
the dependence of the observed fluxes on P and P is 
similar . In fact , before dithering, the pseudo-luminosity 
law of IGonthier et all (j2002f ) scales as P _1 P 1 / 3 , which 
is close to the flux depend ence oc p -° s p 0A in the case 
of the geometric model of IGonthier et all (I2004D derived 
in section 14.4.21 In both cases, the modelled luminosity 
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contour lines are not parallel to the death line and the 
pulsars do not uniformly dim as they approach it. 

4.6. Distribution of Braking Indexes 

As mentioned in section [3. 3. 21 we have assumed in our 
simulations the braking index of each synthetic pulsar to 
be equal to 3, the value for pure magnetic dipole brak¬ 
ing in a vacuum, despite the measured values (c.f. Table 
11 being systematically lower. To demonstrate that our 
results are in fact robust to this assumption, we have 
repeated our “optimal” simulations, but with a normal 
distribution of braking indexes with mean 2.4 and stan¬ 
dard deviation 0.6, the sample estimates for the mea¬ 
sured values, and with an uniform distribution in the 
range [1.4,3.0], bracketing the measured values. The cor¬ 
responding P — P diagrams for typical Monte Carlo re¬ 
alizations of the models were found to be qualitatively 
very similar to the one obtained with all braking indexes 
set to 3 (Figure [7]). This is somewhat surprising, since 
the flow of the pulsars in the P — P diagram depends 
on the braking indexes, and reflects the dominant role of 
the luminosity dependence on P and P in shaping the 
diagram observed in our simulations. 

4.7. Birthrate and Number of Pulsars 

Table [8] lists the derived Galactic pulsar birthrate and 
number of pulsars (potentially observable and total) av¬ 
eraged over the 50 MC realizations of our optimal pop¬ 
ulation synthesis model, calculated assuming the TM98 
beaming model. The formal uncertainties provided are 
the corresponding standard deviations and do not include 
the uncertainty on the model parameters, which is diffi¬ 
cult to reliably quantify. Varying the model parameters 
about the optimal values suggests that the true uncer¬ 
tainties exceed the formal ones by a factor ~ 5. 

The pulsar birthrate estimate, 2.8 ±0.1 pulsars 
per century, is high er than the one derived by 
lYranesevic et al.l (|2004h using a more model-independent 
pulsar current analysis of the PM data (also based on 
the NE2001 and TM98 mode l s) of 1 .9 ± 0.4 pulsars per 
century. As lYranesevic et al.l (j2004l ) explain, their value 
is more strictly interpreted as a lower limit on the true 
value, so that our results are consistent. 

Between 13% and 25% of core-collapse supernovae 
(Type lb, Ic, and II) are expected to leave a black- 
hole remnant instead of a neutron star, depending on 
the lower mass limit for fall-back black hole form atio n, 
for main-sequence stars of solar metallicity dHeger et al.1 
l2003f) . Under the hypothesis that every core-collapse su¬ 
pernova produces either a pulsar or a black hole, the 
expected Galactic core-collapse supern ova r ate is then 
~ 3.2 — 3.7 supernovae per century. iDraeicevich et alJ 
(1199911 reviews estimates based on extragalactic data. 
The values range from ~ 1 to ~ 4 per century, assum¬ 
ing that about 85% of supernovae are core-collapse (e.g., 
iTammann et al.l 1 1994T) . Thus, the data are consistent 
with all neutron stars, except for a small fraction affected 
by accretion from a binary companion (X-ray binaries) 
or endowed with an extraordinarily strong magnetic field 
(magnetars), being born as radio pulsars, although not 
necessarily beamed toward us. Unfortunately, the re¬ 
maining uncertainties prevent us from making a conclu¬ 
sive statement. 


5. CONCLUSION 

Motivated by recent advances in pulsar astronomy (a 
large, homogeneous sample of detections by the Parkes 
and Swinburne Multibeam pulsar surveys; an updated 
model of the Galactic free electron density, NE2001; and 
new astrometric measurements), we have revisited the 
problem of the birth and evolution of isolated radio pul¬ 
sars. 

We started by estimating the pulsar birth velocity 
di stribution dire c tly from pro per motion measurements 
of iBrisken et al.l ( 20021 12003 1 in section [2j A single¬ 
component Maxwellian distribution appears to be inade¬ 
quate due to the detection of a few very high-velocity 
objects. However, we do not find evidence for mul¬ 
timodality of the velocity distribution, as alternative 
single-parameter models with heavier tails accommo¬ 
date the observations equally well as a two-component 
Maxwellian. The exact shape of the velocity distribution 
is not well constrained. We adopted a model in which 
the absolute one-dimensional birth velocity components 
are exponentially distributed and with three-dimensional 
mean (u 3 p) = 380+60 km s -1 . 

In section[3j we used this velocity distribution as input 
to a more general Monte Carlo-based population synthe¬ 
sis. We described parametric prescriptions for the birth 
properties (location, velocity, spin period, magnetic field, 
and radio luminosity) of the pulsars and their time evo¬ 
lution (spatial, rotational, and radio emission shut-off). 
We then generated synthetic pulsar populations on which 
we performed simulations of the PM and SM surveys. By 
comparing the observed samples in the simulations with 
the real detections, we determined “optimal” model pa¬ 
rameters. 

The Galactocentric radial distribution of pulsar for¬ 
mation a ppears consistent w i th the functional form pro¬ 
posed by lYusifov fc Kiicilkl tl2004l f. which incorporates 
a deficit in surfa ce density near t he G alactic Cen¬ 
ter. Although the lYusifov fc Kiiciikl (|2004f) distribution 
was derived for the present-day distribution of evolved 
pulsars, our simulations suggest that younger pulsars 
are preferentially detec ted, a bias not treated in the 
lYusifov fc Kiicilk |2004f ). and that it may in fact be a 
better approximation to the birth distribution. This dis¬ 
tribution is qualitatively consistent (peaking midway be¬ 
tween the Galactic Center and the Sun) with that of mas¬ 
sive main-sequence stars, the purported pulsar progeni¬ 
tors, although there is a difference of ~ 1.5 kpc between 
the radii of peak density, possibly due to remaining im¬ 
perfections in the free electron density model. Spiral arm 
structure is required to reproduce the spatial distribution 
of pulsars. We have modelled this structure using fixed 
spirals, but there are apparent deficiencies. Proper mod¬ 
elling of the angular motion of the spirals with respect 
to the Sun may be needed for more accurate results. 

Pulsar radio luminosities independent of the period 
and period derivative can safely be ruled out. They 
lead to the detection of too many old synthetic pul¬ 
sars, as indicated by the exceedingly large observed scale 
height and a clear pile-up of detections on the death 
line. A model in which the radio (pseudo-)luminosity is, 
before dithering, proportional to p- 15 p°- 5 (the square 
root of the spin-down luminosity) favors the detection 
of younger objects and appears consistent with the ob- 
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TABLE 8 

Galactic Birthrate and Number of Pulsars 


Birthrate 
(psr century -1 ) 

Potentially Observable Pulsars a 

All Radio-Loud Pulsars 

2.8 ±0.1 

120 , 000 ± 20, 000 

1,200, 000 ±300, 000 


Note. — The values given are averages over the 50 Monte Carlo realizations of 
the optimal model. The uncertainties are the corresponding standard deviations 
and do not include the uncertainty on the model parameters, which is difficult 
to reliably quantify. Varying the model parameters about the optimal values 
suggests that the true uncertainties exceed the formal ones by a factor ~ 5. 
a Radio-loud and beamed toward us. 


servations. In this model, the undithered luminosity is 
proportional to the voltage drop available for particle 
acceleration in magnetosphere models, which may be re¬ 
lated to the physics of the radio emission mechanism. 
Also, the modelled fluxes are independent of the viewing 
geometry, provided that the pulsar beam crosses the line 
of sight, which supports the hypothesis of sharp-edged 
beams with random component locations. 

We do not find evidence for significant torque decay 
(due to magnetic field decay or otherwise) over the life¬ 
time of the pulsars as ra dio sources (~ 1 0 0 My r ). Th e 
conflicting conclusion of iGonthier et al.l (120021 . I2004D . 
who found evidence for magnetic field decay on time 
scales < 5 Myr, is most likely due to different assumed 
(effective) dependences of the radio luminosity on the pe¬ 
riod and period derivative of the pulsar, which resulted in 
a pile-up of synthetic pulsars on the death line without 
simulated field decay. Our preferred luminosity model 
neatly avoids such a pile-up, as the undithered luminosity 
contours are parallel to the death line, so that the pulsars 
uniformly dim as they approach it. While it is unclear 
whether we have identified the correct radio pulsar lumi¬ 
nosity law, this argues that conclusions regarding torque 
decay in isolated neutron stars based on simulation stud¬ 
ies are strongly dependent on assumptions regarding the 
luminosities of the pulsars. We have demonstrated that 
it is possible to avoid the conclusion of magnetic field 
decay using a very simple and natural luminosity model. 

Our method does not precisely constrain the distri¬ 
bution of birth spin periods of the pulsars, because the 
period of a pulsar is asymptotically independent of its 
initial value. However, in order to avoid an excess of 
simulated detections with periods < 100 ms, many pul¬ 
sars must be formed with greater periods. We found that 
a Gaussian distribution with mean 300 ms and standard 
deviation 150 ms accommodates the observations well. 
This suggests that characteristic ages, which assume that 
the initial period is negligible compared to the observed 
one, are generally poor age indicators, at least for young 
pulsars. 

We estimate the Galactic pulsar birthrate to be ~ 2.8 
pulsars per century. After accounting for a fraction of 
core-collapse supernovae forming black holes, this value 
is consistent with all neutron stars, except for a small 
fraction in X-ray binaries and magnetars, being born as 
radio pulsars, although not neces sarily b eamed toward 
us. U sing the beaming model of iTauris fe Manchester! 
dl998T) . the Galaxy is estimated to contain ~ 120,000 
potentially observable ordinary pulsars. 


An important point to emphasize is that, as is also the 
case with previous work, we have merely provided ev¬ 
idence for consistency of the observational data with a 
particular scenario. The value of our study lies in the 
fact that our proposed model has, arguably, minimal 
complexity and does not involve any controversial com¬ 
ponent. We have also implemented the recent advances 
in pulsar astronomy in a generally self-consistent man¬ 
ner. That very similar studies differing principally by 
assumptions about poorly constrainted aspects - such as 
the radio luminosity of pulsars - led to conflicting con¬ 
clusions - notably regarding torque decay - is strongly 
suggestive that pulsar population simulations require fur¬ 
ther independent input before they can be used to draw 
definitive conclusions. Perhaps the single most impor¬ 
tant breakthrough would be an independently verified 
luminosity model, as it is what ultimately determines 
the detectability of pulsars. Until then, Ockham’s razor 
should be applied when considering the results of such 
studies. Meanwhile, the pulsar data set is poised to con¬ 
tinue to be enhanced by the on-going Arecibo L-Band 
Feed Array (ALFA) pulsar survey, which is expected to 
detect as m any as 1 000 ordinary pulsars a nd 50 millisec¬ 
ond o nes (jFaucher-Giguere fc Ka spi 2004: G ordes et al.l 
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APPENDIX 

APPENDIX 1: Smin CALCULATION 

We describ e the calculation of the minimum flux theoretically detectable from a radio pulsar ( S m i n ) using the 
iDewev et al.l (S1 9841 ) formula. First, we calculate the observed pulse width, 


vv e = 


(Al) 


where W = DC x P is th e intrinsic pulse widt h, t sarnp is the sampling interval, and T sca tt is the pulse broadening due 
to interstellar scattering dLorimer et al JIT 9931) . DC is the duty cycle, i.e. the fraction of the time that the pulsar’s 
flux is above 50% of its maximum. We take DC = 5%, a typical measured value (see, e.g., Lorimer et al. 2005), for all 
pulsars. We do not model the observed scatter about this value or its potential correlation with the spin period. We 
do, however, model the related period dependence of the beaming fraction the spin period as explained in section [3.41 
Tsamp is the effective sampling interval, taking into account details of the system hardware such as anti-aliasing filters. 
We take T samp = 1.5 t samp . The proportionality coefficient adopted is somewhat arbitrary, but chosen to be of order 
unity. The simulations are insensitive to the exact value of this parameter. N c h is the number of channels across the 
receiver bandwidth and DMq is the dispersion measure at which the smearing of the pulsar in one channel is equal to 
tsamp ■ The latter is given by 


DM 0 = 


Ncht sampV 

8299A u c h 


(A2) 


where DMq is in units of pc cm 3 , t sa mp is in s, and v (the observing frequency) and Av (the receiver bandwidth) 
are in MHz f Johnstonl!l990D . Finally, 


, _ £ l3&{T rec + T sky ) I W e 

'min — Obeam ,. rrrr ?—* : \ t t r 

GyJN p Avt int V P - We 


(A3) 


where, as in section 13.5.21 T rec is the receiver temperature on cold sky, T sky is the sky background temperature, G 
is the antenna gain, N p is the number of polarizations, Av is the receiver bandwidth, t lnt is the integration time, P 
is the pulse period, W e is the effective pulse width, a is the signal-to-noise threshold, and f3 is a constant accounting 
for various sy stem losses. The sky background temperature at 408 MHz is taken from an elec tronic version of th e 
lHaslam et al.l ( 1981 ) map. The value is scaled to 1.4 GHz assuming a spectral index a,b g = —2.8 (lLawson et alill987 ). 
The flux degradation factor due to a pulsar lying away from the center of a telescope beam, assuming a Gaussian 
power pattern, is 

$beam = ^ x p[ {t/ fbeam) \ (A4) 


where fbeam = FWHM /(2v / ln2) (jLorimer et al.lll993 1. In the simulations, we assume that each synthetic pulsar in 
the sky area covered by the PM and SM surveys lies at a random position within the half-power cone of a telescope 
beam. The gain G is t aken to be the ave r age of over individual beams of the Parkes Multibeam receiver, the values of 
which are reported by iManchester et all (120011 ). The values of the parameters for the PM and SM surveys are given 
in Table [5j 


APPENDIX 2: COMPUTATION OF THE BIRTH RADIAL DISTRIBUTION OF PULSARS FROM AN EVOLVED MODEL 

In this appendix, we describe a method to estimate the birth Galactic radial distribution of pulsars corresponding to 
a prescribed mod el fo r t he di stribution of the evolved pulsars. We then apply it to the distribution of evolved pulsars 
derived bv lYusifov fc Kiiclikl (|2004l) . 

We first note that the relationship between the birth and evolved spatial distributions of pulsars depends on the 
details of their evolution (not only spatial, but also rotational, as this determines the pulsars that survive as radio 
sources in the evolved population) and that one can be estimated from the other only after an evolution model has 
been fixed. We thus perform this calculation only after having identified an “optimal” population synthesis model (see 
section ET51) . 

Let p{r ev ) be a prescribed probability density function (PDF) for the evolved radial coordinate r ev of a pulsar, 
i.e. suppose that a pulsar has probability p(r ev )dr ev of having an evolved radial coordinate in the range [r ev ,r ev + 
dr ev ). Given an evolution model, we can compute an “evolution kernel” K(r ev ,rbirth ) by Monte Carlo using the 
computer code described in section [3] Specifically, for equidistant radial distances rbirth in the range [0, 50) kpc 
(fbirth = 0.25 + 0.5* kpc; i = 0,.., 99), we generate MC pulsars with birth radial coordinate equal to rbirth and evolve 
them in time, until 10,000 evolved radio-loud pulsars are reached. We then construct a histogram (with bins centered 
on the selected rbirth values and a bin width of 0.5 kpc) to estimate the distribution of evolved radial distances r ev . 
Interpolating between the discrete values of rbirth and r ev , we obtain a kernel K(r ev , rbirth) defined on the square 
[0, 50) x [0, 50) kpc describing the conditional PDF that a pulsar born with radial distance rbirth has radial distance 
r ev at the time of observation. By the law of total probability, the PDF for the birth radial distance, p{rbirth ), is 
related to p(r ev ) by 

POO 

pij’ev') = / p{j'birth)K{j'ev')^'birth) ( ^'birth' (Bl) 
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Given a birth radial distribution, it is thus straightforward compute the corresponding evolved distribution. Here, 
the evolved distribution is prescribed. We can estimate the birth distribution by choosing a parametric functional 
form, p(rbirth', 9), depending on the vector parameter 9, and carrying out a least squares minimization between the 
prescribed evolved distribution, p(r ev ;presc), and the one corresponding to the trial birth distribution, p(r ev ;trial), 
evaluated at the bin centers. 

We considered two functional forms for p{rbirth', 9). First, a functional form 


p(rbirth', YK04) oc m r th ( ? th 1 

V ft© + fti 


a 

exp 



Tbirth Re \ 

Rq + f?i / 


identical to the one form proposed by lYusifov fe Kiiciikl (l2004f ) for the evolved distribution 10 , 
two-parameter “displaced Gaussian” 


(B2) 

Second, a simpler 


p{r birth', Gauss ) oc exp 


ipbirth Rpeak ) 

^rurth 


(B3) 


truncated to Hirth > 0. Both functional forms resulted in equal least sums of squares, but the fit using the 
p(pbirth', y Tv 04) function was found to be of lesser numerical stability, leading to degenerate and unexpectedly large val¬ 
ues of the parameters in order to reproduce the requisite relatively narrow concentration of pulsars at radius ~ 7 kpc. 
We thus adopted the displaced Gaussian fit ( R pea k = 7.04 and cr rbirth = 1-83). This distribution must not be confused 
with the Gaussian radial distribution discussed in section 14.1.11 In the latter case, the radial distribution in surface 
density is Gaussian and peaks at the Galactic Center. Here, it is the radial PDF p{rbirth'. Gauss) that is Gaussian 
and the peak of the distribution is displaced from the GC. Figure [5] shows the best fit birth radial PDF. To test 
that it satisfactorily reproduces the prescribed evolved distribution, we generated and evolved MC pulsars using our 
optimal population synthesis model and plotted the resulting evolved distribution against the prescribed model. The 
agreement is qualitatively good. 
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